module etienne_gino
use menu_f90
!use S_fitting   ! GINO CONFLICT
use build_lattice   

implicit none
  public
  integer,private,parameter::ndim_old=3
!  integer,private,parameter::sp=kind(1e0)
!  integer,private,parameter::dp=selected_real_kind(2*precision(1e0_sp))
  private graph_data_actlist,graph_data_int
  real(sp) red, green, blue 
  integer, parameter :: nplot=12,n_info=30,N_OBJECTIVE_max=10
  integer, parameter :: nstate=20
  integer :: colour_marker=1001
  integer :: colour_plot=1023
  integer :: colour_layout=1022
  integer :: colour_magnet=1021
  integer :: colour_bend=1020
  integer :: colour_cavity=1013
  integer :: colour_sol=1000
  integer :: colour_wig=998
  integer :: colour_bb=999

  integer :: colour_stable=1019
  integer :: colour_unstable=1018
  integer :: colour_mul(0:3)=(/1014,1017,1016,1015/)
  integer :: colour_random=1013
  integer, pointer :: survey_graphics,xpxgraphics
  integer :: n_turns=1
  integer :: n_ND=3,nfam_mis=3,nfam_scaling=2
  integer, parameter :: nfam_mis_max=5,nfam_scaling_max=2
  integer :: plot(nplot)=0, tplot,tstate,this_plot,form_plot(nplot)=0
  type (GLIMIT) :: plot_size(nplot)
integer :: n_mul=1  ! for plot skipping resonance
!  type (GLIMIT) :: plot_size_sig(nplot)
  character(24) :: plot_sizes_name(nplot)
logical, private :: myorbitunit
integer, parameter :: zturns_nrays = 2
  integer coord(4)
  integer state(nstate)
  logical :: around_fix =.false., store_fix0=.true.,track_a=.false.,track_map=.false.
  logical :: random_colour=.false.
!  real(sp) :: sig_plot(6)=(/1.0,1.0,1.0,1.0,1.0,1.0/)
  integer this_scale
  real :: scale_info(0:n_info)
  real(dp) sig_mis(nfam_mis_max,0:6)
  integer  field_sig_mis(nfam_mis_max,0:6)
  character(1) fam_name(nfam_mis_max)
  character*255 scriptfilename
  character*255 scriptfilenamelist(10)
  character*255 fileflatold
  character*255 fileflatnew
  character*255 trackfile
  character*255 :: lattice_function_file="1_1_1.DAT"
  real(dp) :: y_mis_plot=1.d0
  integer N_OBJECTIVE
  integer  FIELD_OBJECTIVE(N_OBJECTIVE_max)
  real(dp),target ::  OBJECTIVE(N_OBJECTIVE_max)
  real(dp),allocatable ::  zturns(:,:,:)
  real(dp) ztotal_tunes(3,zturns_nrays)
  real(dp) scaling_fact(nfam_scaling_max)
  integer  field_scaling_fact(nfam_mis_max)
  character(2) scaling_fact_name(nfam_mis_max)
  real(dp) previous_ray(6)
  real :: dxmag=0.005
 real :: dxframe=0.05
  type(worm) vers
  LOGICAL :: OSCIL_=.TRUE.,orbit_track=.false.,new_flat=.true.,dirdotted=.false.
  logical :: authorwarn=.true.
type(fibre), pointer :: p_gino_1,p_gino_2
integer :: select_gino_12=1
real(dp) ent_gino(3,3)
real(dp) :: a_gino(3)=0.d0, ang_gino(3)=0.d0
!type (GLIMIT)  :: viewport = GLIMIT(0.0,1.0,0.0,1.0)
type(fibre), pointer :: pbpm
integer::  ibpm=0, wgbpm=0,wibpm=0,wnbpm=0, plane_bpm=1
character(nlp) :: nbpm="#####"
type (GLIMIT3)  :: world = GLIMIT3(-30.0,30.0,-30.0,30.0,-30.0,30.0)
   complex(dp)  am_res(3,zturns_nrays)
   real(dp) am_amp(3,zturns_nrays),rsca(zturns_nrays)
   integer but1,but2
LOGICAL ::  KEEP_COLOUR=.FALSE.  
  type graph_data
   real(sp) s(2*ndim_old)
   real(dp) r(2*ndim_old)
   real(dp) a(6,6)
   LOGICAL  IN(ndim_old),do_a_for_dyn
   integer nd
   integer, pointer :: tc  ! tracking node
   integer nturns,plane,turn_tracked
   integer colour,widget(3)
   logical stable,interrupt,aperture
  end type graph_data

  INTERFACE assignment (=)
     MODULE PROCEDURE graph_data_actlist
     MODULE PROCEDURE graph_data_int
  end  INTERFACE

  type beamlet
   real(dp), pointer :: x(:,:)
   logical, pointer:: s(:)
   integer n
   integer lost
  end type beamlet



contains

 subroutine warn_t(fib,wid)
 implicit none
 type(fibre),pointer :: fib
 character*255 line
 integer wid,mid

mid=fib%mag%p%nst
mid=mid+4

!write(6,*) "associated(fib%parent_layout%t) ", associated(fib%parent_layout%t)

if(associated(fib%parent_layout%t))then

if(mod(mid,2)==0) then

  mid=(mid-4)/2+2
 line=' '
 write(line,'(a30,1x,i4)')  "The mid point is at position ",fib%t1%pos+mid  
 else
  line='No mid point, odd number of steps: click HELP'
endif

else
 line="No Node Layout: no internal tracking"
endif  

 call gmSetTextSetting(wid,line)

 end  subroutine warn_t
 
 subroutine fill_t(fib,w)
 implicit none
 integer w(:)
 type (graph_data)  :: ray
 type(fibre),pointer :: fib
 
  call gmSetValueSetting(w(1),fib%mag%p%nst+4 )
if(associated(fib%parent_layout%t))then
  call gmSetValueSetting(w(2),fib%t1%pos )
  call gmSetValueSetting(w(3),fib%t2%pos )
  call gmSetValueSetting(w(4),fib%t1%pos)  !+fib%mag%p%nst/2+2 )
  call gmSetWidgetProp(w(5) ,gmMinimum =fib%t1%pos, gmMaximum= fib%t2%pos+1  )
  call gmSetWidgetProp(w(5) ,gmTicks1 =fib%mag%p%nst+4, gmTicks2= fib%mag%p%nst+4 )
  call gmSetGaugeValue(w(5) ,fib%t1%pos)  !+fib%mag%p%nst/2+2 )
!  call gmSetWidgetProp(w(5) ,gmIcon=w(6))
else
  call gmSetValueSetting(w(2),0 )
  call gmSetValueSetting(w(3),0 )
  call gmSetValueSetting(w(4),0 )
  call gmSetWidgetProp(w(5) ,gmMinimum =0, gmMaximum= 1 )
  call gmSetGaugeValue(w(5) ,0 )
!  call gmSetWidgetProp(w(5) ,gmIcon=w(8))
endif
!  call gmDrawWidget(w(5))

end subroutine fill_t

 subroutine fill_ray(ray)
 implicit none
  type (graph_data)  :: ray
  real(dp) x(4),b(2,2),bi(2,2),det
  integer i,j

  ! write(6,*) ray%a(1:4,1:4)
   b(1,1)=sqrt(ray%a(1,1)**2+ray%a(1,2)**2)
   b(1,2)=sqrt(ray%a(1,3)**2+ray%a(1,4)**2)
   b(2,2)=sqrt(ray%a(3,3)**2+ray%a(3,4)**2)
   b(2,1)=sqrt(ray%a(3,1)**2+ray%a(3,2)**2)
  !  write(6,*) b

       det=b(1,1)* b(2,2)-b(1,2)*b(2,1)
   bi(1,1)=b(2,2)/det
   bi(2,2)=b(1,1)/det
   bi(1,2)=-b(1,2)/det
   bi(2,1)=-b(2,1)/det
    write(6,*) bi

  
  x=0.0_dp
  x(1)=bi(1,1)*ray%r(1)+bi(1,2)*ray%r(3)
  x(3)=bi(2,1)*ray%r(1)+bi(2,2)*ray%r(3)
  
 !    write(6,*)x
 ray%r(1:4)=0.0_dp
  
  do i=1,4
  do j=1,4
    ray%r(i)=ray%a(i,j)*x(j)+ray%r(i)
  enddo
  enddo
  
!      write(6,*)ray%r(1:4)

 end  subroutine fill_ray


logical function gEnqKeyState_et(key)
implicit none

integer key
if(gEnqKeyState(key)==1) then
 gEnqKeyState_et=.true.
else
 gEnqKeyState_et=.false.
endif

return
end function gEnqKeyState_et

logical function gmEnqToggleSwitch_et(key)
implicit none

integer key
if(gmEnqToggleSwitch(key)==1) then
 gmEnqToggleSwitch_et=.true.
else
 gmEnqToggleSwitch_et=.false.
endif

return
end function gmEnqToggleSwitch_et


subroutine set_gino_warnings
implicit none
 call gSetMaxErrorLimit(-1)   ! no limits on errors
 call gSetErrorMode(galloff)
end subroutine set_gino_warnings

subroutine read_scales(w)

implicit none
integer w
select case(this_scale)
case(0:n_info)
scale_info(this_scale)= gmEnqValueSetting(w )
end select 

end subroutine read_scales

subroutine read_sig_mis
implicit none
integer i,j
do i=1,nfam_mis
do j=0,6
sig_mis(i,j)= gmEnqValueSetting(field_sig_mis(i,j) )
enddo
enddo
end subroutine read_sig_mis

subroutine set_sig_mis

implicit none
integer i,j
do i=1,nfam_mis
do j=0,6
 call  gmSetValueSetting(field_sig_mis(i,j),sig_mis(i,j) )
enddo
enddo

end subroutine set_sig_mis

subroutine read_objective
implicit none
integer i
do i=1,n_objective
 OBJECTIVE(i)= gmEnqValueSetting(FIELD_OBJECTIVE(i) )
enddo
end subroutine read_objective

subroutine set_objective
implicit none
integer i
do i=1,n_objective
 call  gmSetValueSetting(FIELD_OBJECTIVE(i),OBJECTIVE(i) )
enddo

end subroutine set_objective


subroutine print_scales(w)

implicit none
integer w
select case(this_scale)
case(0:n_info)
 call gmSetValueSetting(w,scale_info(this_scale) )
end select 

end subroutine print_scales


subroutine Enq_coordinate_fields

implicit none
 if(this_plot/=0) then
  call gmEnqGraphicsFrameBounds(plot(this_plot),plot_size(this_plot))
!  if(this_plot<=3) then
!   plot_size_sig(this_plot)%xmax=plot_size(this_plot)%xmax/sig_plot(this_plot*2-1)
!   plot_size_sig(this_plot)%xmin=plot_size(this_plot)%xmin/sig_plot(this_plot*2-1)
!   plot_size_sig(this_plot)%ymax=plot_size(this_plot)%ymax/sig_plot(this_plot*2)
!   plot_size_sig(this_plot)%ymin=plot_size(this_plot)%ymin/sig_plot(this_plot*2)
!   else
   plot_size(this_plot)%xmax=plot_size(this_plot)%xmax 
   plot_size(this_plot)%xmin=plot_size(this_plot)%xmin 
   plot_size(this_plot)%ymax=plot_size(this_plot)%ymax 
   plot_size(this_plot)%ymin=plot_size(this_plot)%ymin 
   
!  endif
 endif
end subroutine Enq_coordinate_fields

subroutine set_coordinate_fields

implicit none
if(this_plot/=0) then
  call gmSetValueSetting(coord(1), plot_size(this_plot)%xmin) 
  call gmSetValueSetting(coord(2), plot_size(this_plot)%xmax) 
  call gmSetValueSetting(coord(3), plot_size(this_plot)%ymin) 
  call gmSetValueSetting(coord(4), plot_size(this_plot)%ymax) 
endif
end subroutine set_coordinate_fields

subroutine set_colours

implicit none
      call gDefineRGB(gbackground,0.0, 0.0, 0.0) 
      call gDefineRGB(colour_plot,1.0, 1.0, 1.0) 
      call gDefineRGB(colour_layout,1.0, 0.0, 0.0) 
      call gDefineRGB(colour_magnet,0.0, 1.0, 0.0)  
      call gDefineRGB(colour_bend,0.0, 0.0, 1.0)  
      call gDefineRGB(colour_cavity,0.7, 0.7, 0.0)  
      call gDefineRGB(colour_sol,0.7, 0.0, 0.7)
      call gDefineRGB(colour_wig,1.0, 0.5, 0.0)  
      call gDefineRGB(colour_bb,0.5, 0.1, 0.9)  
      call gDefineRGB(colour_marker,0.0, 0.7, 0.7)  
	  call gDefineRGB(colour_unstable,1.0,0.0,0.0)
	  call gDefineRGB(colour_stable,0.0,1.0,0.0)
	  call gDefineRGB(colour_mul(1),1.0,0.0,0.0)
	  call gDefineRGB(colour_mul(2),0.0,1.0,0.0)
	  call gDefineRGB(colour_mul(3),0.0,0.0,1.0)
	  call gDefineRGB(colour_mul(0),0.7,0.0,0.7)
	  
end subroutine set_colours

recursive subroutine read_pref(file)

!use build_lattice
implicit none
integer mf,i,nf,res,mf0
character(*),optional :: file
character*255 line,filename

call kanalnummer(mf)
if(present(file)) then
 open(unit=mf,file=file,status='old',err=1000)
 filename=file
else
 open(unit=mf,file='pref.txt',status='old',err=1000)
 filename='pref.txt'
endif
read(mf,'(a255)') line
close(mf)

res=index(line,"=")
if(res/=0) then
 filename=line(res+1:255)
 call context(filename)
  write(6,*) filename
 call read_pref(file=filename)
 return
endif

call kanalnummer(mf)
open(unit=mf,file=filename,status='old',err=1000)

read(mf,*,end=100) tplot

do i=1,tplot
! if(i<=3) then
!  read(mf,*) plot_size_sig(i)%xmin,plot_size_sig(i)%xmax &
!  ,plot_size_sig(i)%ymin,plot_size_sig(i)%ymax &
!  ,sig_plot(i*2-1) ,sig_plot(i*2) 
!   plot_size(i)%xmax=plot_size_sig(i)%xmax*sig_plot(i*2-1)
!   plot_size(i)%xmin=plot_size_sig(i)%xmin*sig_plot(i*2-1)
!   plot_size(i)%ymax=plot_size_sig(i)%ymax*sig_plot(i*2)
!   plot_size(i)%ymin=plot_size_sig(i)%ymin*sig_plot(i*2)
! else
  read(mf,*) plot_size(i)%xmin,plot_size(i)%xmax,plot_size(i)%ymin,plot_size(i)%ymax
! endif
enddo

read(mf,*) dxmag, n_nd

READ(MF,*) n_turns


scale_info=100.d0
scale_info(0)=1.d0
 write(6,*) tplot, " frame sizes "
read(mf,'(a120)') trackfile
call context(trackfile)
read(mf,'(a120)') fileflatold
call context(fileflatold)
read(mf,'(a120)') fileflatnew
call context(fileflatnew)


read(mf,*) nf
write(6,*) nf, "script files"
if(nf>10) then
 nf=10
 write(6,*) "Only 10 script files allowed "
endif
do i=1,10
 scriptfilenamelist(i)=' '
enddo
do i=1,nf
 read(mf,*) scriptfilenamelist(i)
 call context(scriptfilenamelist(i))
enddo
scriptfilename=scriptfilenamelist(1)

close(mf)

 ent_gino=0.d0
 a_gino=0.d0 
do i=1,3
 ent_gino(i,i)=1.d0
enddo

tplot=8
allocate(zturns(6,0:n_turns,zturns_nrays))
zturns=zero
ANG_gino=0.d0

return

100 continue

write(6,*) "Empty Preference file "
pause

1000 continue
call kanalnummer(mf0,file="pref.txt")
write(6,*) "No Preference file: defaults used"

tplot=7
do i=1,tplot
plot_size%xmin=-.01
plot_size%xmax=.01
plot_size%ymin=-.001
plot_size%ymax=.001
! if(i<=3) then
!sig_plot(i*2-1)=1.0
!sig_plot(i*2)=1.0
!   plot_size(i)%xmax=plot_size_sig(i)%xmax*sig_plot(i*2-1)
!   plot_size(i)%xmin=plot_size_sig(i)%xmin*sig_plot(i*2-1)
!   plot_size(i)%ymax=plot_size_sig(i)%ymax*sig_plot(i*2)
!   plot_size(i)%ymin=plot_size_sig(i)%ymin*sig_plot(i*2)
!endif
enddo
dxmag=0.01 
n_nd=1
n_turns=1
trackfile="TRACKING.TXT"
fileflatold="flat_old.txt"
fileflatnew="flat_new.txt"
call context(trackfile)
call context(fileflatold)
call context(fileflatnew)
nf=1
do i=1,10
 scriptfilenamelist(i)=' '
enddo
 scriptfilenamelist(1)='script.txt'
call context(scriptfilenamelist(1))


allocate(zturns(6,0:n_turns,zturns_nrays))
zturns=zero

write(mf0,*) tplot

do i=1,tplot
  write(line,*) plot_size(i)%xmin,plot_size(i)%xmax,plot_size(i)%ymin,plot_size(i)%ymax
  write(mf0,'(a255)') line
! endif
enddo
write(mf0,*) dxmag, n_nd

write(mf0,*) n_turns


write(mf0,'(a120)') trackfile
write(mf0,'(a120)')  fileflatold
write(mf0,'(a120)')  fileflatnew


write(mf0,*)  nf
 

do i=1,nf
 write(mf0,*) scriptfilenamelist(i)
enddo


close(mf0)







end subroutine read_pref


subroutine clear_plots

implicit none
integer i
 do i=1,nplot
  if(plot(i)/=0) then 
   call gmClearGraphicsFrame(plot(i))
  endif
 enddo
end subroutine clear_plots

subroutine set_plots_size

implicit none
integer i
 do i=1,nplot
  if(plot(i)/=0) then 
!   if(this_plot<=3) then
!   if(i<=3) then
!   plot_size(this_plot)%xmax=plot_size_sig(this_plot)%xmax*sig_plot(this_plot*2-1)
!   plot_size(this_plot)%xmin=plot_size_sig(this_plot)%xmin*sig_plot(this_plot*2-1)
!   plot_size(this_plot)%ymax=plot_size_sig(this_plot)%ymax*sig_plot(this_plot*2)
!   plot_size(this_plot)%ymin=plot_size_sig(this_plot)%ymin*sig_plot(this_plot*2)
!   plot_size(i)%xmax=plot_size_sig(i)%xmax*sig_plot(i*2-1)
!   plot_size(i)%xmin=plot_size_sig(i)%xmin*sig_plot(i*2-1)
!   plot_size(i)%ymax=plot_size_sig(i)%ymax*sig_plot(i*2)
!   plot_size(i)%ymin=plot_size_sig(i)%ymin*sig_plot(i*2)
!  endif

   call gmSetGraphicsFrameBounds(plot(i),plot_size(i))
  endif
 enddo

 call  gSetViewport3D(world,plot_size(8))
end subroutine set_plots_size

subroutine set_plots_size_special(xm,invert)

implicit none
integer i
real(dp) xm(6)
logical invert
 do i=1,3
   plot_size(i)%xmax=xm(2*i-1)
   plot_size(i)%xmin=-xm(2*i-1)
   plot_size(i)%ymax=xm(2*i)
   plot_size(i)%ymin=-xm(2*i)
   call gmSetGraphicsFrameBounds(plot(i),plot_size(i))
 enddo
if(invert) then
 i=3 
   plot_size(i)%ymax=xm(2*i-1)
   plot_size(i)%ymin=-xm(2*i-1)
   plot_size(i)%xmax=xm(2*i)
   plot_size(i)%xmin=-xm(2*i)
   call gmSetGraphicsFrameBounds(plot(i),plot_size(i))
endif
end subroutine set_plots_size_special

subroutine graph_data_actlist(g,a)

implicit none 
type(graph_data), intent(inout)::g
type (GACTION), intent(in)  :: a
 select case(g%plane)
  case(1:ndim_old) 
  case default
   g%plane=1
  end select
g%s(2*g%plane-1) = a%x
g%s(2*g%plane  ) = a%y
g%r(2*g%plane-1) = a%x
g%r(2*g%plane  ) = a%y

end subroutine graph_data_actlist

subroutine graph_data_int(g,i)
implicit none 
type(graph_data), intent(inout)::g
integer, intent(in)  :: i
 g%s=0.0_dp
 !g%r=0.0_dp
 G%IN=.FALSE.
 g%nturns=i
 g%nd=N_ND
 g%plane=1
 g%stable=.true.
 g%aperture=.true.
 g%interrupt=.false.
 g%colour=1023  ! Gino's last colour user defined
 g%widget=0  ! Usually graphics identifier
end subroutine graph_data_int


subroutine egm_scale_draw_e_beam(b,nd,FIX)
implicit none 
integer nd,i
real x,y
real(dp), optional :: fix(6)
type(e_beam) b     
real f(6)
integer colour
 
f=0.0

if(present(fix)) f=fix
 

!if(nd==1) then
!OSCIL_=.NOT.OSCIL_
!  IF(OSCIL_) THEN
   colour=colour_stable 
!  ELSE
!   colour=gyellow
!  ENDIF  

!endif
 
call gmActivateGraphicsFrame( plot(nd))

!call gmClearGraphicsFrame(plot(nd))

do i=1,b%n

if(.not.b%z(i)%u) then
 
  y=b%z(i)%x(2*nd) - f(2*nd)
  x=b%z(i)%x(2*nd-1) -f(2*nd-1)
 
 call egm_scale_draw(plot(nd),x,y,colour)
endif
if(mod(i,1000)==0) call gFlushGraphics

enddo
 call gFlushGraphics

end subroutine egm_scale_draw_e_beam

subroutine egm_scale_plot_e_beam(b,nd,FIX)
implicit none 
integer nd,i
real x,y,xm,ym
real(dp), optional :: fix(6)
type(e_beam) b     
real f(6)
integer colour
 
f=0.0
xm=0
ym=0
if(present(fix)) f=fix
 


call gmActivateGraphicsFrame( plot(nd))

!call gmClearGraphicsFrame(plot(nd))

do i=1,b%n

if(.not.b%z(i)%u) then
 
  y=abs(b%z(i)%x(2*nd) - f(2*nd))
  x=abs(b%z(i)%x(2*nd-1) -f(2*nd-1))
 
 if(x>xm) xm=x
 if(y>ym) ym=y
endif


enddo

   ym=1.1*ym
   xm=1.1*xm

   plot_size(nd)%xmax=xm
   plot_size(nd)%xmin=-xm
   plot_size(nd)%ymax=ym
   plot_size(nd)%ymin=-ym

if(this_plot==nd) then
  call gmSetValueSetting(coord(1), plot_size(nd)%xmin) 
  call gmSetValueSetting(coord(2), plot_size(nd)%xmax) 
  call gmSetValueSetting(coord(3), plot_size(nd)%ymin) 
  call gmSetValueSetting(coord(4), plot_size(nd)%ymax) 
endif
      call gmSetGraphicsFrameBounds(plot(nd),plot_size(nd))
      call gmDrawWindow(plot(nd) )


 





end subroutine egm_scale_plot_e_beam


subroutine egm_scale_draw_sbeam(b,nd,invert)
implicit none 
integer nd,i
real x,y
type(beamlet) b
logical invert     
integer colour,is

if(nd==1) then
OSCIL_=.NOT.OSCIL_
  IF(OSCIL_) THEN
   colour=colour_stable 
  ELSE
   colour=gyellow
  ENDIF  

endif
 
call gmActivateGraphicsFrame( plot(nd))

call gmClearGraphicsFrame(plot(nd))

is=1
if(b%n>1000) is=b%n/1000

do i=1,b%n,is

if(b%s(i)) then
 if(invert) then
  x=b%x(i,2)
  y=b%x(i,1)
 else
  x=b%x(i,1)
  y=b%x(i,2)
 endif
 call egm_scale_draw(plot(nd),x,y,colour)
!endif
endif
if(mod(i,1000)==0) call gFlushGraphics

enddo
 call gFlushGraphics

end subroutine egm_scale_draw_sbeam


subroutine egm_scale_draw(widget,x,y,colour)

implicit none 
type (GLIMIT)  :: limits
type (GPIXFRAME) :: bounds
real, intent(in):: x,y
integer, intent(in):: widget,colour
integer ix,iy
logical draw
type (GPIXEL) pix
type (GPOINT)  point
      call gmActivateGraphicsFrame( widget)
      DRAW=.TRUE.
      call gmEnqGraphicsFrameBounds(widget, limits) 
      call gmEnqGraphicsPixelBounds(widget, bounds) 
      call gmEnqPixelPos(widget,x,y,pix)
      if(Pix%ix<bounds%xmin.or.Pix%ix>bounds%xmax) draw=.false.
      if(Pix%iy<bounds%ymin.or.Pix%iy>bounds%ymax) draw=.false.
  if(draw) call gDrawPixel(Pix%IX, PIX%IY, colour)

end subroutine egm_scale_draw

!$$$$$$$$$$$$$$$$$$  DAVID BRADLY'S ROUTINES $$$$$$$$$$$$$$$$$$$$$$$$$
subroutine gmEnqPixelPos_et(widget,x,y,pixel)

implicit none
integer, intent(in) :: widget
real, intent(in) :: x,y
type (GPIXEL), intent(out) :: pixel

type (GLIMIT)  :: limits
type (GPIXFRAME) :: bounds
real rx,ry,xi,yi

call gmEnqGraphicsFrameBounds(widget, limits) 
call gmEnqGraphicsPixelBounds(widget, bounds) 
rx=real(bounds%xmax-bounds%xmin)/(limits%xmax-limits%xmin)
ry=-real(bounds%ymax-bounds%ymin)/(limits%ymax-limits%ymin)
xi= (x-limits%xmin)*rx+bounds%xmin
yi= (y-limits%ymin)*ry+bounds%ymax
pixel%ix=int(xi)
pixel%iy=int(yi)
end subroutine gmEnqPixelPos_et

subroutine gmEnqPosofPixel_et(widget,ix,iy,point)

implicit none
integer, intent(in) :: widget
integer, intent(in) :: ix,iy
type (GPOINT), intent(out) :: point

type (GLIMIT)  :: limits
type (GPIXFRAME) :: bounds
real rx,ry,xi,yi

call gmEnqGraphicsFrameBounds(widget, limits) 
call gmEnqGraphicsPixelBounds(widget, bounds) 
rx=(limits%xmax-limits%xmin)/real(bounds%xmax-bounds%xmin)
ry=(limits%ymax-limits%ymin)/(-real(bounds%ymax-bounds%ymin))
xi= (ix-bounds%xmin)*rx+limits%xmin
yi= (iy-bounds%ymin)*ry+limits%ymax
point%x=xi
point%y=yi
return
end subroutine gmEnqPosofPixel_et

!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$   End $$$$$$$$$$$$$$$$$$$$$

subroutine tune_zturns
implicit none
integer i,j,nh,k
type(real_8) y(6)
type(damap) id
type(normalform) normal
real(dp) a(6,6),ztunes0(3,zturns_nrays),ztunes(3,zturns_nrays), &
zdtunes0(3,zturns_nrays),zdtunes(3,zturns_nrays),sca,x(6)
integer, allocatable :: je(:)

my_fix(5)=my_delta
 CALL FIND_ORBIT(my_ering,MY_FIX,MY_START,my_estate,1.D-5)
!write(6,*) my_fix

sca=one/zturns_nrays

do k=zturns_nrays,1,-1
rsca(zturns_nrays-k+1)=k*sca
enddo

call init(my_estate,1,0)
call alloc(id)
call alloc(y)
call alloc(normal)
nh=c_%nd2
if(c_%ndpt/=0) nh=nh-2
id=1
y=my_fix+id

call track(my_ering,y,MY_START,my_estate)

normal=y
!write(6,*) normal%tune

a=zero

do i=1,6
a(i,i)=one
enddo 
allocate(je(c_%nd2))

id=normal%a_t**(-1)

do i=1,nh
do j=1,nh
je(j)=1
a(i,j)=id%v(i).sub.je
je(j)=0
enddo
enddo

do k=2,zturns_nrays

  zturns(1:nh,0,k)= (zturns(1:nh,0,1)-my_fix(1:nh))*rsca(k)+my_fix(1:nh)
  do j=nh+1,6
    zturns(j,0,k)=zturns(j,0,1)
  enddo
  x=zturns(1:nh,0,k)
 do i=1,n_turns
  call track(my_ering,x,MY_START,my_estate)
  zturns(1:nh,i,k)=x
 enddo
enddo



do k=1,zturns_nrays 
 do i=0,n_turns
  zturns(:,i,k)=matmul(a,(zturns(:,i,k)-my_fix))
!  write(16,*) zturns(1,i),zturns(2,i)
 enddo
enddo

ztotal_tunes=0.d0

do k=1,zturns_nrays
do i=0,n_turns
  do j=1,nh/2
   if(i==0) then 
     ztunes0(j,k)=atan2(zturns(2*J-1,i,k), zturns(2*J,i,k))
     if(ztunes0(j,k)<0.0_dp) ztunes0(j,k)=ztunes0(j,k)+twopi
   else
     ztunes(j,k)=atan2(zturns(2*J-1,i,k), zturns(2*J,i,k))
     if(ztunes(j,k)<0.0_dp) ztunes(j,k)=ztunes(j,k)+twopi
     zdtunes(j,k)=(ztunes(j,k)-ztunes0(j,k))/twopi
     ztunes0(j,k)=ztunes(j,k)
     if(zdtunes(j,k)<0.0_dp) zdtunes(j,k)=zdtunes(j,k)+1.0_dp
     ztotal_tunes(j,k)=ztotal_tunes(j,k)+zdtunes(j,k)
   endif
  enddo
enddo
enddo
ztotal_tunes=ztotal_tunes/n_turns
write(6,*) ztotal_tunes

call kill(id)
call kill(y)
call kill(normal)
deallocate(je)

end subroutine tune_zturns

subroutine am_res_zturns
implicit none
integer i,j,nh,k
complex(dp) ex,exn
real(dp) pow

nh=c_%nd2
if(c_%ndpt/=0) nh=nh-2
pow=0
do i=1,nh/2
 pow=pow+abs(mres(i))
enddo
pow=pow-2
!k=1
am_res=0.d0

do k=1,zturns_nrays

ex=zero
do i=1,nh/2
 ex=ex+mres(i)*twopi*ztotal_tunes(i,k)
enddo
ex=exp(-i_*ex)

exn=one
do i=0,n_turns
 do j=1,nh/2
  am_res(j,k)=am_res(j,k)+ exn*(zturns(2*J-1,i,k)**2+ zturns(2*J,i,k)**2) !(zturns(2*J-1,i,k)-i_*zturns(2*J,i,k))   !
 enddo
  exn=exn*ex
enddo
!(1d0+cos(step*(mf-maxn2)))
 do j=1,nh/2
  am_res(j,k)=am_res(j,k)/(n_turns+1)/(zturns(2*J-1,0,k)**2+ zturns(2*J,0,k)**2)
enddo
 do j=1,nh/2
   am_amp(j,k)=abs(am_res(j,k))
 enddo
 
enddo  ! k

write(6,*) " amplitudes ", (mres(j),j=1,nh/2)
do k=1,zturns_nrays
 write(6,*) (am_amp(j,k),j=1,nh/2)
enddo
if(pow/=0) then
 write(6,*) " Scaled amplitudes "  
 do k=1,zturns_nrays
 write(6,*)  (am_amp(j,k)/rsca(k)**pow,j=1,nh/2) 
 enddo
endif

write(6,*)( log(am_amp(j,1)/am_amp(j,2))/log(rsca(1)/rsca(2)),j=1,nh/2)


end subroutine am_res_zturns


subroutine e_gm_draw_grid(colour,col_centre,nx,ny,widget)

implicit none


integer i,col

integer, intent(in) :: colour,col_centre
integer, intent(in) ::widget,nx,ny
type(glimit) line
real(sp) dx,dy,x,y
    call gEnqLineColour(col) 

    call gSetLineColour(colour)  
      call gmActivateGraphicsFrame(widget )
    call gmEnqGraphicsFrameBounds(widget,line)
    dx=(line%xmax-line%xmin)/nx
    dy=(line%ymax-line%ymin)/ny

    do i=0,nx
     x=	i*dx+line%xmin
	 y= line%ymin
     call gMoveTo2D(x,y)  
	 y= line%ymax
	  if(i==nx/2) call gSetLineColour(col_centre)
      call gDrawLineTo2D(x,y) 
	  if(i==nx/2) call gSetLineColour(colour)
    enddo
    do i=0,ny
     y=	i*dy+line%ymin
	 x= line%xmin
     call gMoveTo2D(x,y)  
	 x= line%xmax
	  if(i==ny/2) call gSetLineColour(col_centre)
      call gDrawLineTo2D(x,y) 
	  if(i==ny/2) call gSetLineColour(colour)
    enddo
      call gFlushGraphics
    call gSetLineColour(col)  


end subroutine e_gm_draw_grid

!### Routines used in drawing a moving rectangle while resizing####
subroutine e_gDraw_4_LinesTo2D(widget,line,style)

implicit none
integer, intent(in) ::widget,style
type(glimit), intent(in) :: line
type(glimit) lim
integer p_style



       call gEnqBrokenLine(p_style) 

       call gSetBrokenLine(style) 

       call e_g_order(line,lim)
       call gmActivateGraphicsFrame(widget )
       call gMoveTo2D(lim%xmin,lim%ymax)
       call gDrawLineTo2D(lim%xmax,lim%ymax)
       call gDrawLineTo2D(lim%xmax,lim%ymin)
       call gDrawLineTo2D(lim%xmin,lim%ymin)
       call gDrawLineTo2D(lim%xmin,lim%ymax)
       call gFlushGraphics

       call gSetBrokenLine(p_style) 


end subroutine e_gDraw_4_LinesTo2D

subroutine e_g_order(linein,lim)

implicit none
type(glimit), intent(in) :: linein
type(glimit), intent(out) ::  lim
type(glimit) line
line=linein
                   if(line%xmax>line%xmin) then
                   lim%xmax=line%xmax
                   lim%xmin=line%xmin
                  else
                   lim%xmax=line%xmin
                   lim%xmin=line%xmax
                  endif
                   if(line%ymax>line%ymin) then
                   lim%ymax=line%ymax
                   lim%ymin=line%ymin
                  else
                   lim%ymax=line%ymin
                   lim%ymin=line%ymax
                  endif

end subroutine e_g_order

subroutine e_gscale_g_e_2(g,xg,yg,e,xe,ye)

implicit none
type (GLIMIT)   :: g,e
real xg,yg,xe,ye
xe=xg
ye=yg
return
xe=(e%xmax-e%xmin)/(g%xmax-g%xmin)*(xg-g%xmin) + e%xmin
ye=(e%ymax-e%ymin)/(g%ymax-g%ymin)*(yg-g%ymin) + e%ymin

end subroutine e_gscale_g_e_2

subroutine e_gm_frame_to_pixel_bounds(widget,limits,bounds)

implicit none
integer, intent(in) :: widget
type (GLIMIT), intent(in)  :: limits
type (GPIXFRAME), intent(out)  ::  bounds
type (GLIMIT) lim
type (GPIXEL)  pixel

  call e_g_order(limits,lim)

 call gmEnqPixelPos(widget,lim%xmin,lim%ymax,pixel)
 bounds%xmin=pixel%ix
 bounds%ymin=pixel%iy
 call gmEnqPixelPos(widget,lim%xmax,lim%ymin,pixel)
 bounds%xmax=pixel%ix
 bounds%ymax=pixel%iy

end subroutine e_gm_frame_to_pixel_bounds

!### end ####


!  Gino + PTC
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$   flat survey $$$$$$$$$$$$$$$$$$$$$
!*************************************************************************************
!include "ptc_stuff.f90"
 ! TYPE INTERNAL_STATE
 !    LOGICAL(lp) TOTALPATH,TIME,RADIATION,NOCAVITY,FRINGE,EXACTMIS
 !    LOGICAL(lp) PARA_IN,ONLY_4D,DELTA


subroutine set_ptc_states(my_state)
!!use S_fitting

implicit none
type(internal_state) my_state

 if(my_state%TOTALPATH==1) then ;  call gmSetToggleSwitch(state(1),gon)
 else;                          call gmSetToggleSwitch(state(1),goff); endif
 if(my_state%TIME) then ;  call gmSetToggleSwitch(state(2),gon)
 else;                          call gmSetToggleSwitch(state(2),goff); endif
 if(my_state%RADIATION) then ;  call gmSetToggleSwitch(state(3),gon)
 else;                          call gmSetToggleSwitch(state(3),goff); endif
 if(my_state%NOCAVITY) then ;  call gmSetToggleSwitch(state(4),gon)
 else;                          call gmSetToggleSwitch(state(4),goff); endif
 if(my_state%FRINGE) then ;  call gmSetToggleSwitch(state(5),gon)
 else;                          call gmSetToggleSwitch(state(5),goff); endif
 if(ALWAYS_EXACTMIS) then ;  call gmSetToggleSwitch(state(6),gon)
 else;                          call gmSetToggleSwitch(state(6),goff); endif
 if(my_state%stochastic) then ;  call gmSetToggleSwitch(state(14),gon)
 else;                          call gmSetToggleSwitch(state(14),goff); endif
 if(my_state%ONLY_4D) then ;  call gmSetToggleSwitch(state(7),gon)
 else;                          call gmSetToggleSwitch(state(7),goff); endif
 if(my_state%DELTA) then ;  call gmSetToggleSwitch(state(8),gon)
 else;                          call gmSetToggleSwitch(state(8),goff); endif
 if(around_fix)         then ;  call gmSetToggleSwitch(state(9),gon)
 else;                          call gmSetToggleSwitch(state(9),goff); endif
 if(store_fix0)         then ;  call gmSetToggleSwitch(state(10),gon)
 else;                          call gmSetToggleSwitch(state(10),goff); endif
 if(track_a)            then ;  call gmSetToggleSwitch(state(11),gon)
 else;                          call gmSetToggleSwitch(state(11),goff); endif
 if(c_%APERTURE_FLAG)            then ;  call gmSetToggleSwitch(state(12),gon)
 else;                           call gmSetToggleSwitch(state(12),goff); endif
 if(track_map)            then ;  call gmSetToggleSwitch(state(13),gon)
 else;                          call gmSetToggleSwitch(state(13),goff); endif
 if(random_colour)            then ;  call gmSetToggleSwitch(state(15),gon)
 else;                          call gmSetToggleSwitch(state(15),goff); endif
 if(my_state%spin)            then ;  call gmSetToggleSwitch(state(18),gon)
 else;                          call gmSetToggleSwitch(state(15),goff); endif
 if(new_flat)            then ;  call gmSetToggleSwitch(state(17),gon)
 else;                          call gmSetToggleSwitch(state(17),goff); endif

end subroutine set_ptc_states


subroutine read_ptc_states(my_state)
!!use S_fitting

implicit none
type(internal_state) my_state
integer i

my_state=default0
!around_fix=.false.
if(gmEnqToggleSwitch_et(state(1)))  my_state=my_state+TOTALPATH0
if(gmEnqToggleSwitch_et(state(2)))  my_state=my_state+TIME0
if(gmEnqToggleSwitch_et(state(3)))  my_state=my_state+RADIATION0
if(gmEnqToggleSwitch_et(state(4)))  my_state=my_state+NOCAVITY0
if(gmEnqToggleSwitch_et(state(5)))  my_state=my_state+FRINGE0
if(gmEnqToggleSwitch_et(state(6)))  ALWAYS_EXACTMIS=.true.
if(gmEnqToggleSwitch_et(state(7)))  my_state=my_state+ONLY_4D0
if(gmEnqToggleSwitch_et(state(8)))  my_state=my_state+DELTA0
if(gmEnqToggleSwitch_et(state(14))) my_state=my_state+stochastic0
if(gmEnqToggleSwitch_et(state(18))) my_state=my_state+spin0
if(gmEnqToggleSwitch_et(state(16))) then
 myorbitunit=my_ering%t%ORBIT_LATTICE%ORBIT_USE_ORBIT_UNITS
 orbit_track=.true.
 my_ering%t%ORBIT_LATTICE%ORBIT_USE_ORBIT_UNITS=.false.
 x_orbit_sync(6)=xsm0t
 xsmt=xsm0t
else
 orbit_track=.false.
 if(associated(my_ering%t)) then
  if(associated(my_ering%t%ORBIT_LATTICE)) then
   my_ering%t%ORBIT_LATTICE%ORBIT_USE_ORBIT_UNITS=myorbitunit
  endif
 endif
endif
if(gmEnqToggleSwitch_et(state(17))) then
 new_flat=.true.
else
 new_flat=.false.
endif
call gmSetToggleSwitch(state(9),goff)
!around_fix=.false.
store_fix0=.false.
track_a=.false.
track_map=.false.
c_%APERTURE_FLAG=.false.
if(gmEnqToggleSwitch_et(state(10)))  store_fix0=.true.
if(gmEnqToggleSwitch_et(state(11)))  track_a=.true.
if(gmEnqToggleSwitch_et(state(12)))  c_%APERTURE_FLAG=.true.
if(gmEnqToggleSwitch_et(state(13)))  track_map=.true.

call set_ptc_states(my_state)


end subroutine read_ptc_states

subroutine check_ptc_exceptions(r,MYRING,POS,fix)
!!use S_fitting

implicit none
INTEGER REPORT,I,REPORT_ARRAY(10)
type(graph_data),optional :: r
TYPE(LAYOUT),optional, POINTER::MYRING
INTEGER,optional:: POS
real(dp),optional :: fix(6)
if(present(r)) then
 if(r%interrupt) then
 r%interrupt=.false.
 write(6,*) "interrupted"
 endif
endif

CALL PRODUCE_APERTURE_FLAG(REPORT)
write(6,*) " report ",report
IF(REPORT/=0 ) THEN
 if(present(r))  r%stable=.false.
 CALL RESET_APERTURE_FLAG

! CALL ANALYSE_APERTURE_FLAG(REPORT,REPORT_ARRAY)
 

! WRITE(6,'(A7,1X,10(1X,I4))') "REPORT ",(REPORT_ARRAY(I),I=1,10)

 !WRITE(6,*) "LOST IN MAGNET ",LOST_FIBRE%MAG%NAME
 !CALL FIND_POS(MYRING, LOST_FIBRE ,POS )
 !WRITE(6,*) "This Magnet is magnet number ",pos
 if(present(fix)) fix=zero
ENDIF


if(.not.c_%stable_da) then
write(6,*) "unstable in dabnew"
c_%stable_da=.true.
  if(present(r)) r%stable=.false.
endif

end subroutine check_ptc_exceptions



subroutine e_gm_rescale2d(widget,r)
!!use S_fitting

implicit none
type(layout) r
real(dp) xmax,xmin,zmax,zmin,xave,zave,rd,fac
integer, intent(in)::widget
type(fibre), pointer :: c
integer i
type(glimit) g

nullify(c);

xmax=-10000000.d0
xmin=10000000.d0
zmax=-10000000.d0
zmin=10000000.d0

    C=>r%START
    do i=1,r%n
    if(C%CHART%f%a(1)>xmax) xmax=C%CHART%f%a(1)
    if(C%CHART%f%a(1)<xmin) xmin=C%CHART%f%a(1)
    if(C%CHART%f%a(3)>zmax) zmax=C%CHART%f%a(3)
    if(C%CHART%f%a(3)<zmin) zmin=C%CHART%f%a(3)
    if(C%CHART%f%B(1)>xmax) xmax=C%CHART%f%B(1)
    if(C%CHART%f%B(1)<xmin) xmin=C%CHART%f%B(1)
    if(C%CHART%f%B(3)>zmax) zmax=C%CHART%f%B(3)
    if(C%CHART%f%B(3)<zmin) zmin=C%CHART%f%B(3)
    C=>C%NEXT
    enddo


xave=(xmax+xmin)/2.d0
zave=(zmax+zmin)/2.d0
if(xmax-xmin>zmax-zmin) then
rd=(xmax-xmin)/2
else
rd=(zmax-zmin)/2 
endif

write(6,*) " Dimension of Plot "
write(6,*) xmax,xmin,zmax,zmin,xave,zave
write(6,*) xmax-xmin,zmax-zmin

fac=1.2d0
g%xmin=zave - rd*fac    ! Z
g%ymin=xave - rd*fac 
g%xmax= zave+ rd*fac         !Z
g%ymax= xave+ rd*fac 

call gmSetGraphicsFrameBounds(widget,g)


end subroutine e_gm_rescale2d


subroutine e_gm_rescale2d_universe(widget,mu)
!!use S_fitting

implicit none
type(mad_universe) mu
real(dp) xmax,xmin,zmax,zmin,xave,zave,rd,fac
integer, intent(in)::widget
type(fibre), pointer :: c
integer i,j
type(glimit) g
type(layout), pointer :: r
nullify(c);

xmax=-10000000.d0
xmin=10000000.d0
zmax=-10000000.d0
zmin=10000000.d0
    
    r=>mu%start
    
    do j=1,mu%n
    C=>r%START
    do i=1,r%n
    if(C%CHART%f%a(1)>xmax) xmax=C%CHART%f%a(1)
    if(C%CHART%f%a(1)<xmin) xmin=C%CHART%f%a(1)
    if(C%CHART%f%a(3)>zmax) zmax=C%CHART%f%a(3)
    if(C%CHART%f%a(3)<zmin) zmin=C%CHART%f%a(3)
    if(C%CHART%f%B(1)>xmax) xmax=C%CHART%f%B(1)
    if(C%CHART%f%B(1)<xmin) xmin=C%CHART%f%B(1)
    if(C%CHART%f%B(3)>zmax) zmax=C%CHART%f%B(3)
    if(C%CHART%f%B(3)<zmin) zmin=C%CHART%f%B(3)
    C=>C%NEXT
    enddo
    r=>r%next
    enddo

xave=(xmax+xmin)/2.d0
zave=(zmax+zmin)/2.d0
if(xmax-xmin>zmax-zmin) then
rd=(xmax-xmin)/2
else
rd=(zmax-zmin)/2 
endif

write(6,*) " Dimension of Plot "
write(6,*) xmax,xmin,zmax,zmin,xave,zave
write(6,*) xmax-xmin,zmax-zmin

fac=1.2d0
g%xmin=zave - rd*fac    ! Z
g%ymin=xave - rd*fac 
g%xmax= zave+ rd*fac         !Z
g%ymax= xave+ rd*fac 

call gmSetGraphicsFrameBounds(widget,g)


end subroutine e_gm_rescale2d_universe

subroutine e_gm_draw_layout(r,colour_layout,widget)
!!use S_fitting

implicit none

type(fibre), pointer :: c
integer i,col
type(layout), intent(in) :: r
integer, intent(in) :: colour_layout
integer, intent(in) ::widget
!type(glimit) line
logical draw    !,e_gDrawLineTo2D
!external e_gDrawLineTo2D
nullify(c);

    call gEnqLineColour(col) 

    call gSetLineColour(colour_layout)  
      call gmActivateGraphicsFrame(widget )
!      call gmClearGraphicsFrame(widget )

    C=>r%START
    call gMoveTo2D(real(C%CHART%f%a(3)),real(C%CHART%f%a(1)))  
!    line%xmin=real(C%CHART%f%a(3))
!    line%ymin=real(C%CHART%f%a(1))

    do i=1,r%n
           call gDrawLineTo2D(real(C%CHART%f%a(3)),real(C%CHART%f%a(1))) 
           call gDrawLineTo2D(real(C%CHART%f%b(3)),real(C%CHART%f%b(1))) 
    C=>C%NEXT
    enddo
      call gFlushGraphics
    call gSetLineColour(col)  


end subroutine e_gm_draw_layout

subroutine e_gm_draw_magnets(r,colour_magnet,widget)

!use S_fitting
implicit none
type(layout), intent(in) :: r
type(fibre), pointer :: c
integer i,col
integer, intent(in) :: colour_magnet
integer, intent(in) ::widget
type(glimit) line
logical draw  !,WTF   !,e_gDrawLineTo2D
!external e_gDrawLineTo2D
nullify(c);

! WRITE(6,*) "NEW ?"
!read(5,*) wtf
    call gEnqLineColour(col) 

    call gSetLineColour(colour_magnet)  
      call gmActivateGraphicsFrame(widget )

    C=>R%START
    do i=1,r%n
 !    if(wtf) then
      call e_gm_draw_magnet(widget,c)
  !   else
  !    call e_gm_draw_magnet1(widget,c)
   !  endif
    C=>C%NEXT
    enddo
!WRITE(6,*) C%PREVIOUS%PARENT_LAYOUT%NAME,C%PREVIOUS%MAG%NAME
!WRITE(6,*) C%PARENT_LAYOUT%NAME,C%MAG%NAME
      call gFlushGraphics
    call gSetLineColour(col)  


end subroutine e_gm_draw_magnets

subroutine e_gm_magnet_centre(r)

!use S_fitting
implicit none
type(layout), intent(in) :: r
type(fibre), pointer :: c
integer i
real(dp) lc



a_gino=0.d0

lc=0.d0

    C=>R%START
    do i=1,r%n
      lc=c%mag%p%lc+lc 
       a_gino=C%CHART%f%o*c%mag%p%lc+a_gino
    C=>C%NEXT
    enddo
    
   a_gino=a_gino/lc
   write(6,*) a_gino
end subroutine e_gm_magnet_centre

subroutine e_gm_draw_frame(widget,x,ent,colour)

!use build_lattice
implicit none

integer,optional, intent(in) :: colour
integer, intent(in) ::widget
integer col,i,colm,colb,n,m  
real(dp) x(3),ent(3,3),dx
real  xs(3),xc(3), xa(2),xb(2),phi,dh
type(glimit) line

    call gmEnqGraphicsFrameBounds( widget,line )
    if(line%xmax-line%xmin<line%ymax-line%ymin) then
     dx=(line%xmax-line%xmin)*dxframe
    else
     dx=(line%ymax-line%ymin)*dxframe
    endif

colm=colour_sol  !colour_magnet
if(present(colour)) then
 colm=colour
endif
    call gEnqLineColour(col) 

    call gSetLineColour(colm)  

    call gmActivateGraphicsFrame(widget )

if(present(colour)) then
   call gSetLineColour(colm) 
else
    call gSetLineColour(gorange)  
endif

xs=x(1:3) 
 call project_gino(xs,xa)
call gMoveTo2D(xa(1),xa(2))

   
    xs=dx*ent(3,1:3) + x(1:3) 
 call project_gino(xs,xa)
 call  gDrawArrow2D(xa(1),xa(2),GSOLID)


if(present(colour)) then
   call gSetLineColour(colm) 
else
    call gSetLineColour(gyellow)  ! e_x
endif

 

xs=x(1:3) 
 call project_gino(xs,xa)
call gMoveTo2D(xa(1),xa(2))

   
    xs=dx*ent(1,1:3) + x(1:3) 
 call project_gino(xs,xa)
 call  gDrawArrow2D(xa(1),xa(2),GOPEN)

if(present(colour)) then
   call gSetLineColour(colm) 
else
    call gSetLineColour(ggreen)  ! e_y
endif

xs=x(1:3) 
 call project_gino(xs,xa)
call gMoveTo2D(xa(1),xa(2))

   
    xs=dx*ent(2,1:3) + x(1:3) 
 call project_gino(xs,xa)
 call  gDrawArrow2D(xa(1),xa(2),GOPEN)

    call gSetLineColour(col)  
call gFlushGraphics
end subroutine e_gm_draw_frame

subroutine e_gm_draw_drum(widget,x,z,r,h,colour)

!use build_lattice
implicit none

integer,optional, intent(in) :: colour
integer, intent(in) ::widget
integer col,i,colm,colb,n,m  
real(dp) x,z,r,h,k
real  xs(3),xc(3), xa(2),xb(2),phi,dh

! hiromi

return


colm=colour_sol  !colour_magnet
if(present(colour)) then
 colm=colour
endif
    call gEnqLineColour(col) 

    call gSetLineColour(colm)  

    call gmActivateGraphicsFrame(widget )
    
m=40
n=20
phi=2*pi/n

xc(1)=x
xc(2)=0.0
xc(3)=z

 xs(1)=r+xc(1)
 xs(2)=xc(2)
 xs(3)=xc(3)
call project_gino(xs,xa)

call gMoveTo2D(xa(1),xa(2))

do i=0,n

 xs(1)=r*cos(i*phi)+xc(1)
 xs(2)=0.0+xc(2)
 xs(3)=r*sin(i*phi)+xc(3)

 call project_gino(xs,xa)
 call  gDrawLineTo2D(xa(1),xa(2))

enddo

call gFlushGraphics

xc(1)=x
xc(2)=h
xc(3)=z

 xs(1)=r+xc(1)
 xs(2)=xc(2)
 xs(3)=xc(3)
call project_gino(xs,xa)

call gMoveTo2D(xa(1),xa(2))

do i=0,n

 xs(1)=r*cos(i*phi)+xc(1)
 xs(2)=0.0+xc(2)
 xs(3)=r*sin(i*phi)+xc(3)

 call project_gino(xs,xa)
 call  gDrawLineTo2D(xa(1),xa(2))

enddo

call gFlushGraphics

call gSetBrokenLine(GDOTTED)

xc(1)=x
xc(2)=h
xc(3)=z

 xs(1)=r+xc(1)
 xs(2)=xc(2)
 xs(3)=xc(3)

dh=h/m/n

call project_gino(xs,xa)

call gMoveTo2D(xa(1),xa(2))
phi=phi*1.3
do i=0,n*m

 xs(1)=r*cos(i*phi)+xc(1)
 xs(2)=0.0+xc(2)-i*dh
 xs(3)=r*sin(i*phi)+xc(3)

 call project_gino(xs,xa)
 call  gDrawLineTo2D(xa(1),xa(2))

enddo
call gFlushGraphics

    call gSetLineColour(gred) 
call gSetBrokenLine(GSOLID)
xc(1)=x
xc(2)=kicker_hiromi
xc(3)=z

 xs(1)=r+xc(1)
 xs(2)=xc(2)
 xs(3)=xc(3)
call project_gino(xs,xa)

call gMoveTo2D(xa(1),xa(2))

do i=0,n

 xs(1)=r*cos(i*phi)+xc(1)
 xs(2)=0.0+xc(2)
 xs(3)=r*sin(i*phi)+xc(3)

 call project_gino(xs,xa)
 call  gDrawLineTo2D(xa(1),xa(2))

enddo

call e_gm_draw_frame(widget,x0_hiromi(1:3),ent_hiromi,colour)

    call gSetLineColour(col)  
call gFlushGraphics
if(.false.) then
    call gSetLineColour(gorange)  

xs=x0_hiromi(1:3) 
 call project_gino(xs,xa)
call gMoveTo2D(xa(1),xa(2))
k=r
   
    xs=k*x0_hiromi(4:6) + x0_hiromi(1:3) 
 call project_gino(xs,xa)
 call  gDrawArrow2D(xa(1),xa(2),GSOLID)


    call gSetLineColour(gyellow)  ! e_x
 

xs=x0_hiromi(1:3) 
 call project_gino(xs,xa)
call gMoveTo2D(xa(1),xa(2))
k=r
   
    xs=k*ent_hiromi(1,1:3) + x0_hiromi(1:3) 
 call project_gino(xs,xa)
 call  gDrawArrow2D(xa(1),xa(2),GOPEN)


    call gSetLineColour(ggreen)  ! e_y
xs=x0_hiromi(1:3) 
 call project_gino(xs,xa)
call gMoveTo2D(xa(1),xa(2))
k=r
   
    xs=k*ent_hiromi(2,1:3) + x0_hiromi(1:3) 
 call project_gino(xs,xa)
 call  gDrawArrow2D(xa(1),xa(2),GOPEN)
endif

end subroutine e_gm_draw_drum

subroutine e_gm_track_hiromi(r,widget,colour)
!!use hiromi, only : track_hiromi_in,s_time_hiromi,global_unit_hiromi
integer,optional, intent(in) :: colour
integer, intent(in) ::widget
integer col,i,colm,colb,j
type(graph_data) r
real(dp) x(3),p(3),z0(6),zf(6),t0
type(work) w

z0=r%r
!!call   s_time_hiromi(z0)
 
!!call  global_unit_hiromi(z0)

t0=r%r(6)

write(6,'(a10,6(1x,g12.5))') " initial ",z0

fill_hiromi=.true.
!!call track_hiromi_in(z0,zf,t0)
fill_hiromi=.false.
write(6,'(a10,6(1x,g12.5))') "   final " ,zf

end subroutine e_gm_track_hiromi

subroutine e_gm_draw_hiromi(w,widget,colour)

!use build_lattice
implicit none
integer,optional, intent(in) :: colour
integer, intent(in) ::widget
integer col,i,colm,colb,n,w 
type(glimit) line
real y,dx,x(2),dxd,r(3),dh,rs(4,2),rh(4,2)
type (GHSVSTY) h1,h2
! logical alt
 
!alt=.false.
 
if(w==0) then
 n=n_hiromi
else
 n=n_hiromi_lost
endif
dxd=1.0

colm=colour_plot




if(present(colour)) then
colm=colour
endif
call gEnqHSV(colm,h1)
call gEnqHSV(gbackground,h2)

if(h1%hue==h2%hue.and.h1%sat==h2%sat.and.h1%value==h2%value) return




    call gEnqLineColour(col) 

    call gSetLineColour(colm)  

    call gmActivateGraphicsFrame(widget )

if(w==0) call gSetLineColour(gyellow)     
if(n/=n_hiromi) call gSetLineColour(gred) 
do i=1,n ! =i  n_hiromi
if(w==0) then
 r=z_hiromi0(1:3,i)
else
 r=z_hiromi(1:3,i)
endif
if(i>=above_kicker.and.w/=0.and.(n==n_hiromi)) then
 call gSetLineColour(ggreen)  
endif
 call project_gino(r,x) 
 if(i==1) then
  call gMoveTo2D(x(1),x(2))
 else
  call  gDrawLineTo2D(x(1),x(2))
 endif
 if(mod(i,100)==0) call gFlushGraphics
enddo



    call gSetLineColour(col)  


end subroutine e_gm_draw_hiromi

subroutine e_gm_draw_magnet(widget,c,colour)

!use build_lattice
implicit none
type(fibre), intent(inout):: c
integer,optional, intent(in) :: colour
integer, intent(in) ::widget
integer col,i,colm,colb,n 
type(glimit) line
real y,dx,x(2),dxd,r(3),dh,rs(4,2),rh(4,2)
type (GHSVSTY) h1,h2
 if(.not.c%mag%plot) then
   write(6,*) "c%mag%plot false check old e_gm_draw_magnet1"
   pause 246
   stop
 endif
n=10
dxd=1.0
if(c%dir==-1.and.dirdotted) then
call gSetBrokenLineMode(Gsoft)
call gSetBrokenLine(GLONGCHAINED)
endif
colm=colour_magnet
if(c%mag%p%b0/=0.0_dp) colm=colour_bend
if(c%mag%kind==kind4.or.c%mag%kind==kind21) colm=colour_cavity
if(c%mag%kind==kindwiggler) colm=colour_wig
if(c%mag%kind==kind1) then
 colm=gwhite
 dxd=0.05
endif
if(associated(c%mag%b_sol)) colm=colour_sol
        if(c%mag%p%b0/=0.d0) then
          if(c%mag%bn(1)==0.d0) colm=gcyan 
        endif

    if(c%mag%kind==kind0) then
      colm=colour_marker
!    return  ! drifts and marker
    endif

if(present(colour)) then
colm=colour
endif
call gEnqHSV(colm,h1)
call gEnqHSV(gbackground,h2)

if(h1%hue==h2%hue.and.h1%sat==h2%sat.and.h1%value==h2%value) return


 
    y=y_mis_plot
    call gmEnqGraphicsFrameBounds( widget,line )
    if(line%xmax-line%xmin<line%ymax-line%ymin) then
     dx=(line%xmax-line%xmin)*dxmag  ! *0.012d0
    else
     dx=(line%ymax-line%ymin)*dxmag !*0.012d0
    endif
     dx=dx*dxd

    call gEnqLineColour(col) 

    call gSetLineColour(colm)  

    call gmActivateGraphicsFrame(widget )
    

r=0.5d0*(C%CHART%f%b-C%CHART%f%a)+C%CHART%f%a
 call project_gino(r,x) 

c%i%pos(2)=x(2)
c%i%pos(1)=x(1)

dh=dx/2
    
 r= C%mag%p%f%a -dh*C%mag%p%f%ent(1,1:3)
 r= r + dh*C%mag%p%f%ent(2,1:3)
 call project_gino(r,rs(1,1:2)) 
  r= r + dx*C%mag%p%f%ent(1,1:3)
 call project_gino(r,rs(2,1:2)) 
   r= r - dx*C%mag%p%f%ent(2,1:3)
 call project_gino(r,rs(3,1:2)) 
   r= r - dx*C%mag%p%f%ent(1,1:3)
 call project_gino(r,rs(4,1:2)) 

 
 

 r= C%mag%p%f%b -dh*C%mag%p%f%exi(1,1:3)
 r= r + dh*C%mag%p%f%exi(2,1:3)
 call project_gino(r,rh(1,1:2)) 
  r= r + dx*C%mag%p%f%exi(1,1:3)
 call project_gino(r,rh(2,1:2)) 
   r= r - dx*C%mag%p%f%exi(2,1:3)
 call project_gino(r,rh(3,1:2)) 
   r= r - dx*C%mag%p%f%exi(1,1:3)
 call project_gino(r,rh(4,1:2)) 

call gMoveTo2D(rs(1,1),rs(1,2))
do i=2,4
 call  gDrawLineTo2D(rs(i,1),rs(i,2))
enddo
 call  gDrawLineTo2D(rs(1,1),rs(1,2))

call gMoveTo2D(rh(1,1),rh(1,2))
do i=2,4
 call  gDrawLineTo2D(rh(i,1),rh(i,2))
enddo
 call  gDrawLineTo2D(rh(1,1),rh(1,2))

do i=1,4
 call gMoveTo2D(rs(i,1),rs(i,2))
  call  gDrawLineTo2D(rh(i,1),rh(i,2))
enddo

call gSetBrokenLine(GSOLID)
    call gSetLineColour(col)  


end subroutine e_gm_draw_magnet


subroutine e_locate_bpm(c)

!use S_fitting
implicit none
type(fibre), target :: c
integer i,col

!write(6,*) " in e_locate_bpm "
if(.not.allocated(monitors)) return

!write(6,*) " looking in e_locate_bpm "


col=0
    do i=1,size(monitors)
      if(associated(monitors(i)%p,c)) then
       pbpm=>c
       ibpm=i
       col=i
      exit
      endif
    enddo
!write(6,*) col
if(col/=0) then
call gmSetTextSetting(wnbpm ,pbpm%mag%name)
call gmSetValueSetting(wibpm ,ibpm)
endif

end subroutine e_locate_bpm

subroutine set_size_bpm
implicit none
integer i,jm
real max

if(.not.allocated(monitors)) return
if(ibpm==0) return

jm=ibpm

  plot_size(9)%xmin=0 
  plot_size(9)%xmax=m_turn*1.1 
 max=0.0
  do i=1,m_turn
  if(abs(monitors(jm)%r(1,i))>max) max=abs(monitors(jm)%r(1,i)) 
  if(abs(monitors(jm)%r(3,i))>max) max=abs(monitors(jm)%r(3,i)) 
  enddo
 
  plot_size(9)%ymin=-1.1*max 
  plot_size(9)%ymax=1.1*max
  
  call gmSetValueSetting(coord(1), plot_size(9)%xmin) 
  call gmSetValueSetting(coord(2), plot_size(9)%xmax) 
  call gmSetValueSetting(coord(3), plot_size(9)%ymin) 
  call gmSetValueSetting(coord(4), plot_size(9)%ymax) 
  
  call set_plots_size

end subroutine set_size_bpm


subroutine draw_bpm_data
implicit none
integer i,jm,colb,col,dn,dx
real max,r1(2),r0(2),sc
type(glimit) line

if(.not.allocated(monitors)) return
if(ibpm==0) return
max=0.01
    call gmActivateGraphicsFrame(wgbpm )
    call gmEnqGraphicsFrameBounds( wgbpm,line )


    call gEnqLineColour(col) 
if(plane_bpm==1) then    
  colb=colour_stable
  else
   colb=colour_unstable
 endif
     call gSetLineColour(colb)  
     

  r0(1)=1
  r0(2)=monitors(ibpm)%r(plane_bpm,1)
  call gMoveTo2D(r0(1),r0(2))

  do i=2,m_turn
  r0(1)=i
  r0(2)=monitors(ibpm)%r(plane_bpm,i)
	 call  gDrawLineTo2D(r0(1),r0(2))
  enddo

  colb=gwhite
  call gSetLineColour(colb)  

    r0(1)=line%xmin
    r0(2)=0.0
    r1(1)=line%xmax
    r1(2)=0.0
     call gMoveTo2D(r0(1),r0(2))
	 call  gDrawLineTo2D(r1(1),r1(2))
    
     dn=10
     dx=0
     do while(dx<=line%xmax) 
     
     dx=dx+dn
     sc=1.0
     if(mod(dx,50)==0) sc=sc+1.0
     if(mod(dx,100)==0) sc=sc+1.0
     
     r0(1)=dx
    r0(2)= (line%ymax- line%ymin)*max*sc 
    r1(1)=dx
    r1(2)= -(line%ymax- line%ymin)*max*sc
     call gMoveTo2D(r0(1),r0(2))
	 call  gDrawLineTo2D(r1(1),r1(2))
    
    enddo

      call gFlushGraphics

    call gSetLineColour(col)  


end subroutine draw_bpm_data

subroutine draw_bpm_pq 
implicit none
integer i,my_colour
real x,y

if(.not.allocated(monitors)) return
if(ibpm==0) return


previous_ray=0.d0
previous_ray(1:4)=monitors(ibpm)%xf(1:4,1)
write(6,'(a9,4(1x,g21.14))') "   BPM   ", previous_ray(1:4)
write(6,'(a9,4(1x,g21.14))') " Tracked ", monitors(ibpm)%r(1:4,1)
previous_ray(5)=my_delta
   call gmActivateGraphicsFrame(plot(1))
   do i=1,m_turn
    x=monitors(ibpm)%xf(1,i)
    y=monitors(ibpm)%xf(2,i)
    call egm_scale_draw(plot(1),x,y,colour_plot)  
   enddo
   call gFlushGraphics
   call gmActivateGraphicsFrame(plot(2))
   do i=1,m_turn
    x=monitors(ibpm)%xf(3,i)
    y=monitors(ibpm)%xf(4,i)
    call egm_scale_draw(plot(2),x,y,colour_plot)  
   enddo
   call gFlushGraphics

   call gmActivateGraphicsFrame(plot(6))
   do i=1,m_turn
    x=monitors(ibpm)%xn(1,i)
    y=monitors(ibpm)%xn(2,i)
    call egm_scale_draw(plot(6),x,y,colour_plot)  
   enddo
   call gFlushGraphics
   call gmActivateGraphicsFrame(plot(7))
   do i=1,m_turn
    x=monitors(ibpm)%xn(3,i)
    y=monitors(ibpm)%xn(4,i)
    call egm_scale_draw(plot(7),x,y,colour_plot)  
   enddo
   call gFlushGraphics

end subroutine draw_bpm_pq


subroutine draw_bb(c,widget)

!use build_lattice
implicit none
type(integration_node), pointer :: c
integer, intent(in) ::widget
type(glimit) line
real y,dx,r1(2),dxd
integer col
real(dp) o(3),mid(3,3)

dxd=10.0


    y=y_mis_plot
    call gmEnqGraphicsFrameBounds( widget,line )
    if(line%xmax-line%xmin<line%ymax-line%ymin) then
     dx=(line%xmax-line%xmin)*dxmag  ! *0.012d0
    else
     dx=(line%ymax-line%ymin)*dxmag !*0.012d0
    endif
     dx=dx*dxd

    call gEnqLineColour(col) 

    call gSetLineColour(colour_bb)  

    call gmActivateGraphicsFrame(widget )
    

           CALL INVERSE_FIND_PATCH(c%a,c%ent, c%bb%d,c%bb%a,o,mid)


! write(6,*) "dxd"
!read(5,*) dxd
!o(1)=o(1)+dxd
 r1(1)=o(3)-dx*mid(1,3)
 r1(2)=o(1)-dx*mid(1,1)
 call gMoveTo2D(r1(1),r1(2))
 r1(1)=o(3)
 r1(2)=o(1)
  call  gDrawLineTo2D(r1(1),r1(2))

    call gSetLineColour(gred)  
 r1(1)=o(3)+dx*mid(1,3)
 r1(2)=o(1)+dx*mid(1,1)
  call  gDrawLineTo2D(r1(1),r1(2))

 

!      call gFlushGraphics

    call gSetLineColour(col)  


end subroutine draw_bb

subroutine e_gm_draw_magnet_siamese(widget,c,colour,n)

!use build_lattice
implicit none
type(fibre), pointer:: c
integer  colour(2)
integer, intent(in) ::widget
integer col,i,colm,colb,n
type(glimit) line
real y,dx,r1(2),r2(2),s1(2),s2(2)
type(element), pointer :: mag,magn
type(fibre), pointer:: cn
logical doit
integer :: k=0
real(dp) ent(3,3),a(3)
if(mod(n,2)==0) colm=colour(1)
if(mod(n,2)==1) colm=colour(2)
colb=ggreen

mag=>c%mag
magn=>mag
doit=.true.

    call gmEnqGraphicsFrameBounds( widget,line )

    if(line%xmax-line%xmin>line%ymax-line%ymin) then
     dx=(line%xmax-line%xmin)*dxmag   ! *0.012d0
    else
     dx=(line%ymax-line%ymin)*dxmag  !*0.012d0
    endif


do while(doit) 

cn=>magn%parent_fibre

if(n<0) then
 call e_gm_draw_magnet(widget,cn)
 if(associated(magn%siamese_frame)) then
       CALL INVERSE_FIND_PATCH(magn%p%F%a,magn%p%F%ENT, &
            magn%SIAMESE_FRAME%D,magn%SIAMESE_FRAME%ANGLE,a,ent)
!    call gEnqLineColour(col) 
 
 !   call gSetLineColour(colb)  

 !   call gmActivateGraphicsFrame(widget )
!      r1(1)=a(3)-dx*ent(1,3)
!      r1(2)=a(1)-dx*ent(1,1)
!     call gMoveTo2D(r1(1),r1(2))
!      r1(1)=a(3)+dx*ent(1,3) 
!      r1(2)=a(1)+dx*ent(1,1)
!	 call  gDrawLineTo2D(r1(1),r1(2))
!      r1(1)=a(3)-dx*ent(3,3) 
!      r1(2)=a(1)-dx*ent(3,1)
!     call gMoveTo2D(r1(1),r1(2))
!      r1(1)=a(3)+dx*ent(3,3) 
!      r1(2)=a(1)+dx*ent(3,1)
!	 call  gDrawLineTo2D(r1(1),r1(2))
!	  call gFlushGraphics
!    call gSetLineColour(col)  
call e_gm_draw_frame(widget,a,ent)
 endif
else
 call e_gm_draw_magnet(widget,cn,colm)
 if(associated(magn%siamese_frame)) then
       CALL INVERSE_FIND_PATCH(magn%p%F%a,magn%p%F%ENT, &
            magn%SIAMESE_FRAME%D,magn%SIAMESE_FRAME%ANGLE,a,ent)
!    call gEnqLineColour(col) 
 
!    call gSetLineColour(colb)  

 !   call gmActivateGraphicsFrame(widget )
 !     r1(1)=a(3)-dx*ent(1,3)
 !     r1(2)=a(1)-dx*ent(1,1)
 !    call gMoveTo2D(r1(1),r1(2))
 !     r1(1)=a(3)+dx*ent(1,3) 
 !     r1(2)=a(1)+dx*ent(1,1)
!	 call  gDrawLineTo2D(r1(1),r1(2))
 !     r1(1)=a(3)-dx*ent(3,3) 
 !     r1(2)=a(1)-dx*ent(3,1)
 !    call gMoveTo2D(r1(1),r1(2))
 !     r1(1)=a(3)+dx*ent(3,3) 
 !     r1(2)=a(1)+dx*ent(3,1)
!	 call  gDrawLineTo2D(r1(1),r1(2))
!	  call gFlushGraphics
   call e_gm_draw_frame(widget,a,ent)
k=k+1
if(mod(k,2)==0) then
   call e_gm_draw_frame(widget,a,ent,colour=gbackground)
 
else
   call e_gm_draw_frame(widget,a,ent)
endif
!    call gSetLineColour(col)  
 endif
endif

if(.not.associated(magn%siamese)) exit

magn=>magn%siamese

if(associated(magn,mag)) doit=.false.
enddo

      call gFlushGraphics

    call gSetLineColour(col)  


end subroutine e_gm_draw_magnet_siamese

subroutine e_gm_draw_magnet_girder(widget,c,colour,n)

!use build_lattice
implicit none
type(fibre), pointer:: c
integer  colour(2)
integer, intent(in) ::widget
integer col,i,colm,colb,n
type(glimit) line
real y,dx,r1(2) 
type(element), pointer :: mag,magn
type(fibre), pointer:: cn
logical doit
integer :: k=0
!
if(mod(n,2)==0) colm=colour(1)
if(mod(n,2)==1) colm=colour(2)
colb=gorange
mag=>c%mag
magn=>mag
doit=.true.
    call gmEnqGraphicsFrameBounds( widget,line )

    if(line%xmax-line%xmin>line%ymax-line%ymin) then
     dx=(line%xmax-line%xmin)*dxmag   ! *0.012d0
    else
     dx=(line%ymax-line%ymin)*dxmag  !*0.012d0
    endif

do while(doit) 

cn=>magn%parent_fibre


if(n<0) then
 call e_gm_draw_magnet(widget,cn)
 if(associated(magn%girder_frame)) then
call e_gm_draw_frame(widget,magn%girder_frame%b,magn%girder_frame%exi)
!    call gEnqLineColour(col) 
 
!    call gSetLineColour(colb)  

!    call gmActivateGraphicsFrame(widget )
!      r1(1)=magn%girder_frame%b(3)-dx*magn%girder_frame%exi(1,3) 
!      r1(2)=magn%girder_frame%b(1)-dx*magn%girder_frame%exi(1,1) 
!     call gMoveTo2D(r1(1),r1(2))
!      r1(1)=magn%girder_frame%b(3)+dx*magn%girder_frame%exi(1,3) 
!      r1(2)=magn%girder_frame%b(1)+dx*magn%girder_frame%exi(1,1)!
!	 call  gDrawLineTo2D(r1(1),r1(2))
!      r1(1)=magn%girder_frame%b(3)-dx*magn%girder_frame%exi(3,3) 
!      r1(2)=magn%girder_frame%b(1)-dx*magn%girder_frame%exi(3,1)
!     call gMoveTo2D(r1(1),r1(2))
!      r1(1)=magn%girder_frame%b(3)+dx*magn%girder_frame%exi(3,3) 
!      r1(2)=magn%girder_frame%b(1)+dx*magn%girder_frame%exi(3,1)!
!	 call  gDrawLineTo2D(r1(1),r1(2))
!	  call gFlushGraphics
!    call gSetLineColour(col)  
 endif
else
 call e_gm_draw_magnet(widget,cn,colm)
 if(associated(magn%girder_frame)) then
k=k+1
if(mod(k,2)==0) then
call e_gm_draw_frame(widget,magn%girder_frame%b,magn%girder_frame%exi,colour=gbackground)
else
call e_gm_draw_frame(widget,magn%girder_frame%b,magn%girder_frame%exi)
endif
!    call gEnqLineColour(col) 
! 
!    call gSetLineColour(colb)  
!
!    call gmActivateGraphicsFrame(widget )
!      r1(1)=magn%girder_frame%b(3)-dx*magn%girder_frame%exi(1,3) 
!      r1(2)=magn%girder_frame%b(1)-dx*magn%girder_frame%exi(1,1) 
!     call gMoveTo2D(r1(1),r1(2))
!      r1(1)=magn%girder_frame%b(3)+dx*magn%girder_frame%exi(1,3) 
!      r1(2)=magn%girder_frame%b(1)+dx*magn%girder_frame%exi(1,1)
!	 call  gDrawLineTo2D(r1(1),r1(2))
!      r1(1)=magn%girder_frame%b(3)-dx*magn%girder_frame%exi(3,3) 
!      r1(2)=magn%girder_frame%b(1)-dx*magn%girder_frame%exi(3,1)
!     call gMoveTo2D(r1(1),r1(2))
!      r1(1)=magn%girder_frame%b(3)+dx*magn%girder_frame%exi(3,3) 
!      r1(2)=magn%girder_frame%b(1)+dx*magn%girder_frame%exi(3,1)
!	 call  gDrawLineTo2D(r1(1),r1(2))
!	  call gFlushGraphics
!    call gSetLineColour(col)  
 endif
endif

if(.not.associated(magn%girderS)) exit

magn=>magn%girderS

if(associated(magn,mag)) doit=.false.
enddo

      call gFlushGraphics

    call gSetLineColour(col)  


end subroutine e_gm_draw_magnet_girder

subroutine write_girders(fib,imax)
implicit none
type(fibre), target :: fib
type(element), pointer :: e1,e2
integer, optional :: imax
integer i,im
logical print_info,found,ADDIN
type(girder_info), pointer :: inf
real(dp) a(3),ent(3,3),b(3),c(3),d1,d2,exi(3,3),d(3),ang(3),xmis(6)
real(dp) o(3),DELO(3)
im=10000
print_info=.true.
if(present(imax)) im=imax

if(im<0) then
 print_info=.false.
 im=-im
endif


   if(associated(fib%mag%assembly).and.print_info) then  !!!!! fib%mag%assembly
    write(6,*) "Assembly # ",fib%mag%assembly%i
    write(6,*) " # in original assembly ",fib%mag%assembly%mul
    write(6,*) " # added to assembly ",fib%mag%assembly%added
     inf=>fib%mag%assembly%info
     if(associated(inf%mag,fib%mag)) then
      write(6,*) " found immediately "
         found=.true.
     else
      found=.false.
      do while(associated(inf))
        if(associated(inf%mag,fib%mag)) then
         found=.true.
         exit
        endif
      inf=>inf%next
      enddo
      if(found) then
       write(6,*) " found later "
      else
        write(6,*) " must be a drift or marker : not on list "
      endif

     endif
found=.false.
      if(found) then
       write(6,*) inf%mag%name
       write(6,*) inf%d(2),inf%d(11),inf%d(2)+inf%d(11)
       write(6,*) inf%d(3),inf%d(12),inf%d(3)+inf%d(12)
       write(6,*) inf%d(4),inf%d(13),inf%d(4)+inf%d(13)
       FOUND=.FALSE.
       ADDIN=.FALSE.
       CALL FIND_AFFINE_GIRDER(inf%mag%parent_fibre,e1,FOUND)
       IF(FOUND) then
         CALL FIND_FRAME_GIRDER(e1,a,ent,addin)
       endif
!       call FIND_PATCH(inf%a ,inf%ent,fib%mag%p%f%o,fib%mag%p%f%mid,D,ANG)
       call FIND_PATCH(fib%mag%assembly%a,fib%mag%assembly%ent,a ,ent,D,ANG)

       write(6,*) " find patch translation "
       write(6,*) D
       write(6,*) inf%d(11:13)
       write(6,*) " find patch rotation "
       write(6,*) ang
       write(6,*) inf%d(14:16)
!       call FIND_PATCH(a ,ent,fib%mag%assembly%a,fib%mag%assembly%ent,D,ANG)
!       write(6,*) " find patch inverse translation "
!       write(6,*) D
!       write(6,*) " find patch inverse rotation "
!       write(6,*) ang
       
       DELO=fib%mag%P%F%O -fib%mag%assembly%a   ! O"-OMEGA
       
       DELO=MATMUL(fib%mag%assembly%ent,DELO)   !G(O"-OMEGA)

       CALL GEO_ROTA(ENT,DELO,ANG,-1)
       
       DELO=DELO-D-inf%a+fib%mag%assembly%a

       DELO=MATMUL(fib%mag%assembly%ent,DELO)   !G(O"-OMEGA)
       WRITE(6,*) DELO
      
      endif
   endif !if(associated(inf%mag,fib%mag)) then
i=0
   if(associated(fib%mag%girders)) then
   e1=>fib%mag
   e2=> e1%girders
    if(print_info) write(6,*) e1%name,' ',e1%vorname
   do while(.not.associated(e1,e2))
     if(i>im) exit
     i=i+1
    if(print_info) write(6,*) e2%name,' ',e2%vorname
    e2=>e2%girders
   enddo
   endif
! if(.not.(print_info)) THEN
  IF(I>IM) THEN
   WRITE(6,*) " DISASTER IN A GIRDER "
   WRITE(6,*)  e1%name,' ',e1%vorname
   WRITE(6,*)  e2%name,' ',e2%vorname
  ENDIF
! ENDIF
end subroutine write_girders

subroutine write_siamese(fib,imax)
implicit none
type(fibre), target :: fib
type(element), pointer :: e1,e2
integer, optional :: imax
integer i,im
logical print_info
im=10000
print_info=.true.
if(present(imax)) im=imax

if(im<0) then
 print_info=.false.
 im=-im
endif

i=0
   if(associated(fib%mag%siamese)) then

   e1=>fib%mag
   e2=> e1%siamese
    if(print_info) write(6,*) e1%name,' ',e1%vorname
   do while(.not.associated(e1,e2))
     if(i>im) exit
     i=i+1
    if(print_info) write(6,*) e2%name,' ',e2%vorname
    e2=>e2%siamese
   enddo
   endif

!if(.not.(print_info)) THEN
 IF(I>IM) THEN
  WRITE(6,*) " DISASTER IN A SIAMESE "
  WRITE(6,*)  e1%name,' ',e1%vorname
  WRITE(6,*)  e2%name,' ',e2%vorname
 ENDIF
!ENDIF

end subroutine write_siamese

subroutine check_siamese_girder_topology(U,imax)
implicit none
type(mad_universe), target :: U
type(fibre), pointer :: fib
integer i,imax

fib=>u%start%start
do i=1,u%nf

 call write_siamese(fib,imax)
 call write_girders(fib,imax)

if(mod(i,u%nf/10)==0) then
write(6,*) i, "magnets checked "

endif

fib=>fib%n
enddo

write(6,*) "check_siamese_girder_topology done"

end subroutine check_siamese_girder_topology

subroutine e_gm_draw_magnet_3d(c)

!use build_lattice
implicit none
type(fibre), target :: c

integer widget
integer col,i,colm,colb
type(glimit) line
real dx,dy
real xc,yc,zc,dvx,dvy,dvz,r,d
real r1(3),r2(3),r3(3),r4(3)
real r5(3),r6(3),r7(3),r8(3)

widget=plot(8)

colm=colour_magnet
if(c%mag%p%b0/=0.0_dp) colm=colour_bend
if(c%mag%kind==kind4.or.c%mag%kind==kind21) colm=colour_cavity
if(c%mag%kind==kindwiggler) colm=colour_wig
if(associated(c%mag%b_sol)) colm=colour_sol

xc=0.0
yc=0.0
zc=0.0
dvx=0.0
dvy=0.0
dvz=-1.0
r=10.0
d=5.0
dx=1.d0
dy=1.d0

    if(c%mag%kind==kind0.or.c%mag%kind==kind1) return  ! drifts and marker
        if(c%mag%p%b0/=0.d0) then
          if(c%mag%bn(1)==0.d0) colm=gcyan 
        endif

    call gEnqLineColour(col) 

    call gSetLineColour(colm)  
    call gmActivateGraphicsFrame(widget )

 call gDefineSphericalView(xc,yc,zc,r,dvx,dvy,dvz,d)
 call gUpdateView   

 do i=1,3   
  r1(i)=C%mag%p%f%a(i)+dx*C%mag%p%f%ent(1,i)+dy*C%mag%p%f%ent(2,i)
 enddo
 do i=1,3   
  r2(i)=C%mag%p%f%a(i)-dx*C%mag%p%f%ent(1,i)+dy*C%mag%p%f%ent(2,i)
 enddo
 do i=1,3   
  r3(i)=C%mag%p%f%a(i)-dx*C%mag%p%f%ent(1,i)-dy*C%mag%p%f%ent(2,i)
 enddo
 do i=1,3   
  r4(i)=C%mag%p%f%a(i)+dx*C%mag%p%f%ent(1,i)-dy*C%mag%p%f%ent(2,i)
 enddo

 do i=1,3   
  r5(i)=C%mag%p%f%b(i)+dx*C%mag%p%f%exi(1,i)+dy*C%mag%p%f%exi(2,i)
 enddo
 do i=1,3   
  r6(i)=C%mag%p%f%b(i)-dx*C%mag%p%f%exi(1,i)+dy*C%mag%p%f%exi(2,i)
 enddo
 do i=1,3   
  r7(i)=C%mag%p%f%b(i)-dx*C%mag%p%f%exi(1,i)-dy*C%mag%p%f%exi(2,i)
 enddo
 do i=1,3   
  r8(i)=C%mag%p%f%b(i)+dx*C%mag%p%f%exi(1,i)-dy*C%mag%p%f%exi(2,i)
 enddo

   call e_gm_draw_cube_3d(r1,r2,r3,r4,r5,r6,r7,r8)

   call gFlushGraphics

    call gSetLineColour(col)  


end subroutine e_gm_draw_magnet_3d


subroutine e_gm_draw_cube_3d(r1,r2,r3,r4,r5,r6,r7,r8)
implicit none
real r1(3),r2(3),r3(3),r4(3)
real r5(3),r6(3),r7(3),r8(3)

 call et_gMoveTo3D(r1)
 call et_gDrawLineTo3D(r2)
 call et_gDrawLineTo3D(r3)
 call et_gDrawLineTo3D(r4)
 call et_gDrawLineTo3D(r1)
 call et_gDrawLineTo3D(r5)  ! r1 to r5
 call et_gDrawLineTo3D(r6)
 call et_gDrawLineTo3D(r7)
 call et_gDrawLineTo3D(r8)
 call et_gDrawLineTo3D(r5)

 call et_gMoveTo3D(r2)
 call et_gDrawLineTo3D(r6)
 call et_gMoveTo3D(r3)
 call et_gDrawLineTo3D(r7)
 call et_gMoveTo3D(r4)
 call et_gDrawLineTo3D(r8)

end subroutine e_gm_draw_cube_3d


subroutine et_gDrawLineTo3D(r)
implicit none
real r(3),r1(3)

r1=r
 call convert_gino(r1)

call gDrawLineTo3D(r1(1),r1(2),r1(3))

end subroutine et_gDrawLineTo3D

subroutine et_gMoveTo3D(r)
implicit none
real r(3),r1(3)

r1=r
 call convert_gino(r1)

call gMoveTo3D(r1(1),r1(2),r1(3))

end subroutine et_gMoveTo3D


subroutine convert_gino(r1)
implicit none
real r1(3),r(3)
real(dp) a(3),d(3),ent(3,3)


         d=r1-a_gino
  

        ent=ent_gino

        call change_basis(d,global_frame,a,ent)  

        call GEO_ROT(ent,A,ANG_gino,1,ent_gino)

  
r(1)=a(3) 
r(2)=a(1) 
r(3)=a(2) 

r1=r


  


end subroutine convert_gino

subroutine project_gino(r,r1)
implicit none
real r1(2),r(3)
real(dp) a(3),d(3),ent(3,3)
        

         d=r-a_gino
         

      !  ent=ent_gino

      !  call change_basis(d,global_frame,a,ent)  

      !  call GEO_ROT(ent,A,ANG_gino,1,ent_gino)

  
r1(2)=ent_gino(1,1)*d(1)+ent_gino(1,2)*d(2)+ent_gino(1,3)*d(3)
r1(1)=ent_gino(3,1)*d(1)+ent_gino(3,2)*d(2)+ent_gino(3,3)*d(3)

r1(1)=r1(1)+a_gino(3)
r1(2)=r1(2)+a_gino(1)

end subroutine project_gino

subroutine e_gm_draw_muls(r,colour_magnet,nmul,widget)

!use S_fitting
implicit none
type(layout), intent(in) :: r
type(fibre), pointer :: c
integer i,col
integer, intent(in) :: colour_magnet
integer, intent(in) ::widget
type(glimit) line
logical draw    
real bn_max
integer nmul
nullify(c);

    call gEnqLineColour(col) 

    call gSetLineColour(colour_magnet)  
      call gmActivateGraphicsFrame(widget )
if(nmul>0) then
 bn_max=0.0
    C=>R%START
    do i=1,r%n
     if(c%mag%p%nmul>=nmul) then
      if(abs(c%mag%bn(nmul))>bn_max) bn_max = abs(c%mag%bn(nmul))
     endif
    C=>C%NEXT
    enddo
  if(bn_max==0.0) then
    call gSetLineColour(col)  
   return
  endif
    C=>R%START
    do i=1,r%n
      call e_gm_draw_mul(widget,c,bn_max,nmul)
    C=>C%NEXT
    enddo
      call gFlushGraphics
else


 bn_max=0.0
    C=>R%START
    do i=1,r%n
     if(associated(c%mag%b_sol)) then
      if(abs(c%mag%b_sol)>bn_max) bn_max = abs(c%mag%b_sol)
     endif
    C=>C%NEXT
    enddo
  if(bn_max==0.0) then
    call gSetLineColour(col)  
   return
  endif
    C=>R%START
    do i=1,r%n
      call e_gm_draw_mul(widget,c,bn_max,nmul)
    C=>C%NEXT
    enddo
      call gFlushGraphics

endif
    call gSetLineColour(col)  
   
end subroutine e_gm_draw_muls


subroutine e_gm_draw_mul(widget,c,bn_max,nmul,colour)

!use build_lattice
implicit none

type(fibre), intent(inout):: c
integer,optional, intent(in) :: colour
integer, intent(in) ::widget
integer col,i,colm,colb
type(glimit) line
real y,dx,r1(2),r2(2),s1(2),s2(2)
 real dmul
 real bn_max,ndx,dm
 integer dmuln
 integer nmul
 
 
 colm=colour_magnet
if(present(colour)) then
 colm=colour
endif
 if(nmul>4) stop 999
 colm=colour_mul(nmul)

ndx=8.0 
 dmuln=5
dm=sqrt((C%CHART%f%a(1)-C%CHART%f%b(1))**2+(C%CHART%f%a(3)-C%CHART%f%b(3))**2)/2.d0

if(c%mag%kind==kind0.or.c%mag%kind==kind1) return  ! drifts and marker
if(nmul>0) then


if(c%mag%p%nmul<nmul) return
if(c%mag%bn(nmul)==zero) return
    y=y_mis_plot
    call gmEnqGraphicsFrameBounds( widget,line )
    if(line%xmax-line%xmin>line%ymax-line%ymin) then
     dx=(line%xmax-line%xmin)*dxmag  ! *0.012d0
    else
     dx=(line%ymax-line%ymin)*dxmag !*0.012d0
    endif

    call gEnqLineColour(col) 

    call gSetLineColour(colm)  
    call gmActivateGraphicsFrame(widget )
 
 dmul=0.5+float(nmul-1)/dmuln/2.0
    
 r1(2)=dmul*(C%CHART%f%b(1)-C%CHART%f%a(1))+C%CHART%f%a(1)
 r1(1)=dmul*(C%CHART%f%b(3)-C%CHART%f%a(3))+C%CHART%f%a(3)

    
 call gMoveTo2D(r1(1),r1(2))

 r2(1)=  r1(1)+(c%mag%bn(nmul)*dx*ndx/bn_max)*C%mag%p%f%mid(1,3)
 r2(2)= r1(2)+(c%mag%bn(nmul)*dx*ndx/bn_max)*C%mag%p%f%mid(1,1)
 call  gDrawLineTo2D(r2(1),r2(2))



! s1(1)=C%CHART%f%b(3)+y*(C%mag%p%f%b(3)-C%CHART%f%b(3))-dx*C%mag%p%f%exi(1,3)
! s1(2)=C%CHART%f%b(1)+y*(C%mag%p%f%b(1)-C%CHART%f%b(1))-dx*C%mag%p%f%exi(1,1)
! call gMoveTo2D(s1(1),s1(2))

!s2(1)=C%CHART%f%b(3)+y*(C%mag%p%f%b(3)-C%CHART%f%b(3))+dx*C%mag%p%f%exi(1,3)
!s2(2)=C%CHART%f%b(1)+y*(C%mag%p%f%b(1)-C%CHART%f%b(1))+dx*C%mag%p%f%exi(1,1)
! call  gDrawLineTo2D(s2(1),s2(2))

! call gMoveTo2D(r2(1),r2(2))
! call  gDrawLineTo2D(s2(1),s2(2))

! call gMoveTo2D(r1(1),r1(2))
! call  gDrawLineTo2D(s1(1),s1(2))

else
if(.not.associated(c%mag%b_sol)) return

    y=y_mis_plot
    call gmEnqGraphicsFrameBounds( widget,line )
    if(line%xmax-line%xmin>line%ymax-line%ymin) then
     dx=(line%xmax-line%xmin)*dxmag  ! *0.012d0
    else
     dx=(line%ymax-line%ymin)*dxmag !*0.012d0
    endif

    call gEnqLineColour(col) 

    call gSetLineColour(colm)  
    call gmActivateGraphicsFrame(widget )
 
 dmul=0.5+float(nmul-1)/dmuln/2.0
    
 r1(2)=dmul*(C%CHART%f%b(1)-C%CHART%f%a(1))+C%CHART%f%a(1)
 r1(1)=dmul*(C%CHART%f%b(3)-C%CHART%f%a(3))+C%CHART%f%a(3)

    
 call gMoveTo2D(r1(1),r1(2))

 r2(1)=  r1(1)+(c%mag%b_sol*dx*ndx/bn_max)*C%mag%p%f%mid(1,3)
 r2(2)= r1(2)+(c%mag%b_sol*dx*ndx/bn_max)*C%mag%p%f%mid(1,1)
 call  gDrawLineTo2D(r2(1),r2(2))




endif

    call gSetLineColour(col)  


end subroutine e_gm_draw_mul


subroutine e_gm_draw_ray(colour_plot,widget,c,scale,xi,xf)

!use build_lattice
implicit none

type(fibre), intent(in):: c
integer, intent(in) :: colour_plot
integer, intent(in) ::widget
integer col,i,kn
real xi,xf,scale
real r1(2),r2(2),xr,midx,midz,alpha





    call gEnqLineColour(col) 

    call gSetLineColour(colour_plot)  
    call gmActivateGraphicsFrame(widget )
    

 if(c%dir==1) then   
  r1(1)=C%CHART%f%a(3)+scale*xi*C%mag%p%f%ent(1,3)
  r1(2)=C%CHART%f%a(1)+scale*xi*C%mag%p%f%ent(1,1)
  call gMoveTo2D(r1(1),r1(2))

 
  r2(1)=C%CHART%f%b(3)+scale*xf*C%mag%p%f%exi(1,3)
  r2(2)=C%CHART%f%b(1)+scale*xf*C%mag%p%f%exi(1,1)
  call  gDrawLineTo2D(r2(1),r2(2))


 else
  r1(1)=C%CHART%f%b(3)+scale*xi*C%mag%p%f%exi(1,3)
  r1(2)=C%CHART%f%b(1)+scale*xi*C%mag%p%f%exi(1,1)
  call gMoveTo2D(r1(1),r1(2))
 
  r2(1)=C%CHART%f%a(3)+scale*xf*C%mag%p%f%ent(1,3)
  r2(2)=C%CHART%f%a(1)+scale*xf*C%mag%p%f%ent(1,1)
  call  gDrawLineTo2D(r2(1),r2(2))
 endif 

!      call gFlushGraphics

    call gSetLineColour(col)  


end subroutine e_gm_draw_ray

subroutine e_gm_draw_ray_node(colour_plot,widget,c)

!use build_lattice
implicit none

type(three_d_info), intent(inout):: c
integer, intent(in) :: colour_plot
integer, intent(in) ::widget
integer col,i,kn,j
real scale
real r1(2),r2(2),midx,midz,alpha,x(3)
real(dp) ent(3,3),a(3),xr(6),omega0(3),omega(3)




    call gEnqLineColour(col) 

  OSCIL_=.NOT.OSCIL_
  IF(OSCIL_) THEN
    call gSetLineColour(colour_plot)  
  ELSE
    call gSetLineColour(colour_LAYOUT)  
  ENDIF  
    call gmActivateGraphicsFrame(widget )
    


!     omega(1)= a(1)+scale*(xr(1)*ent(1,1)+xr(3)*ent(2,1))
!     omega(2)= a(2)+scale*(xr(1)*ent(1,2)+xr(3)*ent(2,2))
!     omega(3)= a(3)+scale*(xr(1)*ent(1,3)+xr(3)*ent(2,3))
 !    r1(1)=c%r0(3)
 !    r1(2)=c%r0(1)
x=c%r0
 call project_gino(x,r1) 
 !    r2(1)=c%r(3)
 !    r2(2)=c%r(1)
x=c%r
 call project_gino(x,r2) 

     if(abs(r1(1))>1.d6) r1=r2
     call gMoveTo2D(r1(1),r1(2))
     call  gDrawLineTo2D(r2(1),r2(2))


   !   call gFlushGraphics

    call gSetLineColour(col)  


end subroutine e_gm_draw_ray_node



subroutine e_gm_draw_ray_worm(colour_plot,widget,c,scale,omega0)

!use build_lattice
implicit none

type(fibre), intent(in):: c
integer, intent(in) :: colour_plot
integer, intent(in) ::widget
integer col,i,kn,j
real scale
real r1(2),r2(2),midx,midz,alpha
real(dp) ent(3,3),a(3),xr(6),omega0(3),omega(3)




    call gEnqLineColour(col) 

    call gSetLineColour(colour_plot)  
    call gmActivateGraphicsFrame(widget )
    

 DO J=vers%POS(2)+1,vers%POS(3)-1     ! pos(2)+1 to pos(3)-1 inside the magnet
   call gMID(vers,xr,j)   

     call G_FRAME(vers%e,ENT,A,j)

     omega(1)= a(1)+scale*(xr(1)*ent(1,1)+xr(3)*ent(2,1))
     omega(2)= a(2)+scale*(xr(1)*ent(1,2)+xr(3)*ent(2,2))
     omega(3)= a(3)+scale*(xr(1)*ent(1,3)+xr(3)*ent(2,3))
     r1(1)=omega0(3)
     r1(2)=omega0(1)
     r2(1)=omega(3)
     r2(2)=omega(1)
     if(abs(r1(1))>1.d6) r1=r2
     call gMoveTo2D(r1(1),r1(2))
     call  gDrawLineTo2D(r2(1),r2(2))
     omega0=omega
enddo


   !   call gFlushGraphics

    call gSetLineColour(col)  


end subroutine e_gm_draw_ray_worm


subroutine e_gm_draw_fix(r,fix0,colour_plot,widget,scale)

!use build_lattice
implicit none

type(layout), intent(in) :: r
type(fibre), pointer :: c
integer, intent(in) :: colour_plot
integer, intent(in) ::widget
integer col,i
real xi,xf,scale
real r1(2),r2(2)
logical fix0



    call gEnqLineColour(col) 

    call gSetLineColour(colour_plot)  
    call gmActivateGraphicsFrame(widget )
c=>r%start
do i=1,r%n    

if(fix0) then
 xi=c%i%fix0(1)
 xf=c%next%i%fix0(1)
else
 xi=c%i%fix(1)
 xf=c%next%i%fix(1)
endif

 if(c%dir==1) then   
  r1(1)=C%CHART%f%a(3)+scale*xi*C%mag%p%f%ent(1,3)
  r1(2)=C%CHART%f%a(1)+scale*xi*C%mag%p%f%ent(1,1)
  call gMoveTo2D(r1(1),r1(2))
 
  r2(1)=C%CHART%f%b(3)+scale*xf*C%mag%p%f%exi(1,3)
  r2(2)=C%CHART%f%b(1)+scale*xf*C%mag%p%f%exi(1,1)
  call  gDrawLineTo2D(r2(1),r2(2))
 else
  r1(1)=C%CHART%f%b(3)+scale*xi*C%mag%p%f%exi(1,3)
  r1(2)=C%CHART%f%b(1)+scale*xi*C%mag%p%f%exi(1,1)
  call gMoveTo2D(r1(1),r1(2))
 
  r2(1)=C%CHART%f%a(3)+scale*xf*C%mag%p%f%ent(1,3)
  r2(2)=C%CHART%f%a(1)+scale*xf*C%mag%p%f%ent(1,1)
  call  gDrawLineTo2D(r2(1),r2(2))
 endif 
 c=>c%next
enddo   

      call gFlushGraphics

    call gSetLineColour(col)  


end subroutine e_gm_draw_fix

subroutine e_plot_info(r,c,k)

!use build_lattice
implicit none

type(layout), intent(in) :: r
type(fibre), pointer :: c,p
integer, intent(in) :: k
integer mf,i
real(dp) b

write(6,*) " fibre name ,k=n_lat "
write(6,*)  c%mag%name,k
write(6,*)  lattice_function_file

call kanalnummer(mf)
open(unit=mf,file=lattice_function_file)
rewind mf
111  FORMAT(2(1x,E15.8),1x,a16)

 if(.not.associated(r%t)) call make_node_layout(r)
 p=>c
 do i=1,r%n
  b=p%i%beta(k)
  write(mf,111) p%t1%s(1),b,p%mag%name
  write(6,111) p%t1%s(1),b,p%mag%name
  p=>p%next
 enddo
  b=p%i%beta(k)
    write(6,111) p%t1%s(1),b,p%mag%name
write(mf,111) p%previous%t2%s(1),b,p%mag%name

close(mf)

end subroutine e_plot_info

subroutine e_plot_info_monitors(r,k)

!use build_lattice
implicit none

type(layout), intent(in) :: r
type(fibre), pointer :: p
integer, intent(in) :: k
integer mf,i
real(dp) b

write(6,*) " k=n_lat "
write(6,*)  k
write(6,*)  lattice_function_file

call kanalnummer(mf)
open(unit=mf,file=lattice_function_file)
rewind mf
111  FORMAT(2(1x,E15.8),1x,a16)

 do i=1,size(monitors)
  if(monitors(i)%full) then
   p=>monitors(i)%p
   b=p%i%beta(k)
  write(mf,111) p%t1%s(1),b,p%mag%name
  write(6,111) p%t1%s(1),b,p%mag%name
  endif
 enddo
  b=p%i%beta(k)
    write(6,111) p%t1%s(1),b,p%mag%name
write(mf,111) p%previous%t2%s(1),b,p%mag%name

close(mf)

end subroutine e_plot_info_monitors

subroutine e_gm_draw_info(colour_plot,widget)

!use build_lattice
implicit none

type(fibre), pointer :: c
integer, intent(in) :: colour_plot
integer, intent(in) ::widget
integer col,i
real xi,xf,scale
real r1(2),r2(2),r(3)


c=>my_ering%start

do i=1,my_ering%n
    call gEnqLineColour(col) 

    call gSetLineColour(colour_plot)  
    call gmActivateGraphicsFrame(widget )

 scale=scale_info(this_scale)   
 xi=c%i%beta(this_scale)
 xf=xi

 if(associated(c%next)) then
 xf=c%next%i%beta(this_scale)
   if(this_scale<=3) then
    if(xi>xf) xi= 0.0   ! phase advance rules
   endif
 endif


 if(c%dir==1) then 
  r=  C%CHART%f%a+scale*xi*C%mag%p%f%ent(1,1:3)
   call project_gino(r,r1) 

  call gMoveTo2D(r1(1),r1(2))
 
    r=  C%CHART%f%b+scale*xf*C%mag%p%f%exi(1,1:3)
   call project_gino(r,r2) 

 
 else
    r=  C%CHART%f%b+scale*xi*C%mag%p%f%exi(1,1:3)
   call project_gino(r,r1)  

  call gMoveTo2D(r1(1),r1(2))
 
  r=  C%CHART%f%a+scale*xf*C%mag%p%f%ent(1,1:3)
   call project_gino(r,r2) 
 
 endif 
     call  gDrawLineTo2D(r2(1),r2(2))
c=>c%next
enddo
      call gFlushGraphics

    call gSetLineColour(col)  


end subroutine e_gm_draw_info

subroutine e_gm_draw_info1(colour_plot,widget)

!use build_lattice
implicit none

type(fibre), pointer :: c
integer, intent(in) :: colour_plot
integer, intent(in) ::widget
integer col,i
real xi,xf,scale
real r1(2),r2(2)


c=>my_ering%start

do i=1,my_ering%n
    call gEnqLineColour(col) 

    call gSetLineColour(colour_plot)  
    call gmActivateGraphicsFrame(widget )

 scale=scale_info(this_scale)   
 xi=c%i%beta(this_scale)
 xf=xi

 if(associated(c%next)) then
 xf=c%next%i%beta(this_scale)
   if(this_scale<=3) then
    if(xi>xf) xi= 0.0   ! phase advance rules
   endif
 endif



 if(c%dir==1) then   
  r1(1)=C%CHART%f%a(3)+scale*xi*C%mag%p%f%ent(1,3)
  r1(2)=C%CHART%f%a(1)+scale*xi*C%mag%p%f%ent(1,1)
  call gMoveTo2D(r1(1),r1(2))
 

  r2(1)=C%CHART%f%b(3)+scale*xf*C%mag%p%f%exi(1,3)
  r2(2)=C%CHART%f%b(1)+scale*xf*C%mag%p%f%exi(1,1)
 else
  r1(1)=C%CHART%f%b(3)+scale*xi*C%mag%p%f%exi(1,3)
  r1(2)=C%CHART%f%b(1)+scale*xi*C%mag%p%f%exi(1,1)
  call gMoveTo2D(r1(1),r1(2))
 
  r2(1)=C%CHART%f%a(3)+scale*xf*C%mag%p%f%ent(1,3)
  r2(2)=C%CHART%f%a(1)+scale*xf*C%mag%p%f%ent(1,1)
 endif 
     call  gDrawLineTo2D(r2(1),r2(2))
c=>c%next
enddo
      call gFlushGraphics

    call gSetLineColour(col)  


end subroutine e_gm_draw_info1


subroutine e_gm_locate(widget,x,y,name,mag)

!use build_lattice
implicit none

type(fibre), pointer :: c,this
integer, intent(in) ::widget
integer col,i,mag
real x,y,min,d
character(nlp) name
   TYPE(fibre_appearance), POINTER :: ou

min=100000000.0
c=>my_ering%start

do i=1,my_ering%n
d=abs(c%i%pos(2)-y)+abs(c%i%pos(1)-x)
if(d<min) then
  this=>c
  min=d
  name=c%mag%name
  mag=i
endif
c=>c%next
enddo

call e_gm_draw_magnet(widget,this,colour=gred)
call gFlushGraphics

write(6,*) name,mag,c%mag%p%method,c%mag%p%nst
if(this%mag%p%nmul/=0)  then
 do i=1,this%mag%p%nmul
   write(6,*) i,this%mag%an(i),this%mag%bn(i)
 enddo
endif
if(associated(this%mag%b_sol))  then
   write(6,*) "b_sol = ",this%mag%b_sol
endif
if(associated(this%mag%volt))  then
   write(6,*) "volt = ",this%mag%volt
endif

c=>my_ering%start
do i=1,my_ering%n
if(associated(this%mag,c%mag).and.i/=mag) then
Write(6,*) "Also positioned at ",i
endif
c=>c%next
enddo

if(associated(this%mag%doko)) then
ou=>this%mag%doko
 do while(associated(ou))
  write(6,*) "doko ",ou%parent_fibre%pos
  ou=>ou%next
 enddo
endif


end subroutine e_gm_locate

!call e_gm_locate_pointer(plot(this_plot),actlist%x,actlist%y,name,i)

subroutine e_gm_locate_pointer(widget,x,y,name,mag,drawm)

!use build_lattice
implicit none

type(fibre), pointer :: c,this
integer, intent(in) ::widget
integer col,i,mag,n,nt
real x,y,min,d
character(nlp) name
logical, optional :: drawm
logical draw

draw=.true.
if(present(drawm)) draw=drawm


call TIE_MAD_UNIVERSE(M_u)
call Universe_max_n(n)
if(m_t%n/=0) then
call TIE_MAD_UNIVERSE(M_t)
call gUniverse_max_n(m_t,nt)
endif
WRITE(6,*) " UNIVERSE POPULATION ",N,nt
 
min=100000000.0
c=>my_ering%start

do i=1,max(n,nt)
d=abs(c%i%pos(2)-y)+abs(c%i%pos(1)-x)
if(d<min) then
  this=>c
  min=d
  name=c%mag%name
  mag=i
endif
c=>c%n
!c=>c%next
enddo
if(select_gino_12==1) then
 p_gino_1=>this
else
 p_gino_2=>this
endif

if(draw) then
 call e_gm_draw_magnet(widget,this,colour=gred)
 call gFlushGraphics
endif


write(6,*) " Magnet type is ",mytype(this%mag%kind)

write(6,*) "kind,name,pos "

write(6,*)this%mag%kind, name,mag 

write(6,*) "method,nst,exact " 
write(6,*) this%mag%p%method,this%mag%p%nst,this%mag%p%exact

if(this%mag%p%nmul/=0)  then
 do i=1,this%mag%p%nmul
   write(6,*) i,this%mag%an(i),this%mag%bn(i)
 enddo
endif
if(associated(this%mag%b_sol))  then
   write(6,*) "b_sol = ",this%mag%b_sol
endif
if(associated(this%mag%volt))  then
   write(6,*) "volt = ",this%mag%volt
endif


end subroutine e_gm_locate_pointer


subroutine e_gm_locate_pointer_gird(widget,x,y,name,mag,drawm)

!use build_lattice
implicit none

type(layout), pointer :: l
type(fibre), pointer :: c,this
integer, intent(in) ::widget
integer col,i,mag,n,nt
real x,y,min,d,xs(2),r(3)
character(nlp) name
logical, optional :: drawm
logical draw

TYPE(fibre_appearance), POINTER :: ou

draw=.true.
if(present(drawm)) draw=drawm

call TIE_MAD_UNIVERSE(M_u)
call Universe_max_n(n)

if(m_t%n/=0) then
call TIE_MAD_UNIVERSE(M_t)
call gUniverse_max_n(m_t,nt)
endif
WRITE(6,*) " UNIVERSE POPULATION ",N,nt
min=100000000.0


c=>my_ering%start
!if(.not.associated(c%p%mag%p%f))

do i=1,max(n,nt)
 !write(6,*) i,n,c%mag%name
 r=c%mag%p%f%o
 call project_gino(r,xs) 
d=abs(xs(2)-y)+abs(xs(1)-x)
if(d<min) then
  this=>c
  min=d
  name=c%mag%name
  mag=i
endif
c=>c%n
!c=>c%next
enddo

!if(select_gino_12==1) then
! p_gino_1=>this
!else
! p_gino_2=>this
!endif

if(draw) then
 call e_gm_draw_magnet(widget,this,colour=gred)
 call gFlushGraphics
endif
write(6,*) name,mag,c%mag%p%method,c%mag%p%nst
if(this%mag%p%nmul/=0)  then
 do i=1,this%mag%p%nmul
   write(6,*) i,this%mag%an(i),this%mag%bn(i)
 enddo
endif
if(associated(this%mag%b_sol))  then
   write(6,*) "b_sol = ",this%mag%b_sol
endif
if(associated(this%mag%volt))  then
   write(6,*) "volt = ",this%mag%volt
endif

if(associated(this%mag%doko)) then
ou=>this%mag%doko
 do while(associated(ou))
  write(6,*) "doko ",ou%parent_fibre%pos
  ou=>ou%next
 enddo
endif

end subroutine e_gm_locate_pointer_gird


subroutine e_gm_locate_special(widget,x,y,name,mag)

!use build_lattice
implicit none

type(fibre), pointer :: c,this
integer, intent(in) ::widget
integer col,i,mag
real x,y,min,d,d1,d2,d3
character(nlp) name

min=100000000.0
c=>my_ering%start

do i=1,my_ering%n
d=abs(c%i%pos(2)-y)+abs(c%i%pos(1)-x)
if(d<min) then
  this=>c
  min=d
  name=c%mag%name
  mag=i
endif
c=>c%next
enddo

d1=sqrt((c%mag%p%f%a(3)-x)**2+(c%mag%p%f%a(1)-y)**2)
d2=sqrt((c%mag%p%f%o(3)-x)**2+(c%mag%p%f%o(1)-y)**2)
d3=sqrt((c%mag%p%f%b(3)-x)**2+(c%mag%p%f%b(1)-y)**2)

write(6,*) d1,d2,d3




!call e_gm_draw_magnet(widget,this,colour=gred)
call gFlushGraphics

end subroutine e_gm_locate_special



subroutine e_gm_draw_cross(colour_plot,widget,scale,x,y)

!use build_lattice
implicit none

integer, intent(in) :: colour_plot
integer, intent(in) ::widget
type(glimit) line
integer col,i
real x,y,scale,dh,dv
real r1(2),r2(2)


   call gmEnqGraphicsFrameBounds( widget,line )

dh=abs(scale*(line%xmax-line%xmin))
dv=abs(scale*(line%ymax-line%ymin))


    call gEnqLineColour(col) 

    call gSetLineColour(colour_plot)  
    call gmActivateGraphicsFrame(widget )
    
  call gMoveTo2D(x,y-dv)
  call  gDrawLineTo2D(x,y+dv)
  call gMoveTo2D(x-dh,y)
  call  gDrawLineTo2D(x+dh,y)

      call gFlushGraphics

    call gSetLineColour(col)  


end subroutine e_gm_draw_cross

subroutine e_gm_draw_rect(colour_plot,widget,scale,x,y)

!use build_lattice
implicit none

integer, intent(in) :: colour_plot
integer, intent(in) ::widget
type(glimit) line
integer col,i
real x,y,scale,dh,dv
real r1(2),r2(2)
type (GPOINT) :: pts(4) 

   call gmEnqGraphicsFrameBounds( widget,line )

dh=abs(scale*(line%xmax-line%xmin))
dv=abs(scale*(line%ymax-line%ymin))


    call gEnqLineColour(col) 

    call gSetLineColour(colour_plot)  
    call gmActivateGraphicsFrame(widget )

!


   call gMoveTo2D(x-dh,y-dv)
pts(1)=GPOINT(2.0*dh  ,0.0)
pts(2)=GPOINT(0.0     ,2.0*dv)
pts(3)=GPOINT(-2.0*dh ,0.0)
pts(4)=GPOINT(0.0     ,-2.0*dv)
call gDrawPolylineby2D(4,pts);
call gSetLineColour(colour_plot)  

 call gFillPolygonBy2D(gsolid,gcurrent,garea,4,pts)

      call gFlushGraphics

    call gSetLineColour(col)  

end subroutine e_gm_draw_rect


subroutine e_gm_draw_rect_g(colour_plot,widget,x,y,dh,dv,hash)

!use build_lattice
implicit none

integer, intent(in) :: colour_plot
integer, intent(in) ::widget
type(glimit) line
integer col,i
real x,y,scale,dh,dv ,N,dh0,dv0,h,v
real r1(2),r2(2)
type (GPOINT) :: pts(4) 
integer, optional :: hash
N=1000
    call gEnqLineColour(col) 

    call gSetLineColour(colour_plot)  
    call gmActivateGraphicsFrame(widget )

call gmEnqGraphicsFrameBounds( widget,line )

dh0=((line%xmax-line%xmin)/N)
dv0=((line%ymax-line%ymin)/N)
h=abs(line%xmax-line%xmin)
v=abs(line%ymax-line%ymin)



DO I=1,N
 IF((line%xmin+I*DH0<=x-dh).OR. (line%xmin+I*DH0>=x+dh) )then
   call gMoveTo2D(line%xmin+I*DH0,line%ymin)
   call gdrawlineby2D(0.0,v)
  
 endif
 IF((line%ymin+I*Dv0<=y-dv).OR. (line%ymin+I*Dv0>=y+dv) )then
  call gMoveTo2D(line%xmin,line%ymin+I*Dv0 )
  call gdrawlineby2D(h,0.0)
 endif
ENDDO








      call gFlushGraphics

    call gSetLineColour(col)  



end subroutine e_gm_draw_rect_g



! Misalignments
function box_dis(sig)
!USE DFLIB
implicit none 
 REAL(8) ran
 real(8) box_dis,sig
 box_dis=0.d0

ran=bran(ran)
! CALL RANDOM(ran)
box_dis=ran-0.5

box_dis=sqrt(12.d0)*box_dis*sig

end function box_dis


subroutine misalign(r,fam)           ! misaligns families
!use build_lattice
type(layout), intent(inout):: r
integer fam,i,j
type(fibre), pointer:: p
real(dp) mis(6)

!CALL SEED (int(sig_mis(fam,0)))

p=>r%start
do i=1,r%n
if(p%mag%name(1:1)==fam_name(fam)) then
 do j=1,6
  mis(j)=1.d-3*box_dis(sig_mis(fam,j))
 enddo
 p=mis
endif
p=>p%next
enddo

end subroutine misalign

subroutine scale_fam(r,fam)    ! scales some families (bumps here)
!use build_lattice
type(layout), intent(inout):: r
real(dp) v
integer fam,i,j
type(fibre), pointer:: p



p=>r%start
do j=1,r%n
if(p%mag%name(1:2)==scaling_fact_name(fam)) then
 do i=1,p%mag%p%nmul
  v=p%mag%bn(i)*scaling_fact(fam)
  call add(p,i,0,V)
  v=p%mag%an(i)*scaling_fact(fam)
  call add(p,-i,0,V)
 enddo  

endif
p=>p%next
enddo
  if(scaling_fact(fam)==0.d0) then 
   scaling_fact(fam)=1.d0
  else
   scaling_fact(fam)=0.d0
  endif

end subroutine scale_fam





subroutine store_fix(r) ! stores fixed points in fix or fix0
!use build_lattice
integer i
type(layout), intent(inout)::r
type(fibre), pointer:: p
real(dp) x(6)

call move_to(r,p,my_start )
x=my_fix
write(6,*) x
do i=1,r%n
call track(my_ering,x,i,i+1,my_estate)
if(store_fix0) then
 p%i%fix0=x
else
 p%i%fix=x
endif
 p=>p%next
enddo
write(6,*) x
end subroutine store_fix

subroutine auto_scale_info(r,widget,per,k) ! stores fixed points in fix or fix0
!use build_lattice
integer i,widget,k
type(layout), intent(in)::r
type(fibre), pointer:: p
real(sp) min,max,dx,per
type(glimit) line

    call gmEnqGraphicsFrameBounds( widget,line )
    if(line%xmax-line%xmin>line%ymax-line%ymin) then
     dx=(line%ymax-line%ymin)
    else
     dx=(line%xmax-line%xmin)
    endif


min=1.0e10
max=-min


p=>r%start
do i=1,r%n
 if(abs(p%i%beta(k))>max) max=abs(p%i%beta(k))
 if(abs(p%i%beta(k))<min) min=abs(p%i%beta(k))
 p=>p%next
enddo
 scale_info(k)=1.0
 if(max-min/=0.0) then
   scale_info(k)=per*dx/(max)
 endif

end subroutine auto_scale_info

subroutine store_aperture(r) ! stores fixed points in fix or fix0
!use build_lattice
integer i
type(layout), intent(inout)::r
type(fibre), pointer:: p

p=>r%start
do i=1,r%n
!call remove_aperture(p)
if(p%mag%name==name_ap) then
call put_aperture(p,kind_ap,r_ap,x_ap,y_ap)
 Write(6,*) " Aperture put at element ",i," Name = ", p%mag%name
    Write(6,*) 'kind',   p%mag%p%aperture%kind
   Write(6,*) 'dx'  ,  p%mag%p%aperture%dx 
   Write(6,*) 'dy'  ,  p%mag%p%aperture%dy 
   Write(6,*) 'x'  , p%mag%p%aperture%x
   Write(6,*) 'y'  , p%mag%p%aperture%y
    Write(6,*)'r'  , p%mag%p%aperture%r
    Write(6,*)  c_%check_madx_aperture
    Write(6,*)   c_%aperture_flag
endif
p=>p%next
enddo

end subroutine store_aperture

subroutine remove_all_aperture(r) ! stores fixed points in fix or fix0
!use build_lattice
integer i
type(layout), intent(inout)::r
type(fibre), pointer:: p

p=>r%start
do i=1,r%n
call remove_aperture(p)
p=>p%next
enddo

end subroutine remove_all_aperture

subroutine scan_aperture(r,w1,w2) ! stores fixed points in fix or fix0
!use build_lattice
integer i,w1,w2
type(layout), intent(inout)::r
type(fibre), pointer:: p
real x1,x2
x1=-1000
x2=-1000

p=>r%start
do i=1,r%n
if(p%mag%name==name_ap.and.associated(p%mag%p%aperture)) then
 x1=p%mag%p%aperture%x
 x2=p%mag%p%aperture%y
 exit
endif
p=>p%next
enddo

x_ap=x1
y_ap=x2
call gmSetValueSetting(w1,x1)
call gmSetValueSetting(w2,x2)

end subroutine scan_aperture


      integer function frasav(widget,file)
!
! Routine to save contents of graphics frame into an image file
!
! Input arguments:
!     widget : integer : graphics frame widget identifier
!     file   : character : exported image file name
!            The following file extensions are catered for:
!                   bmp, jpeg, png, eps and wmf.
!
! Return Values:
!     integer : 0 = failure
!               1 = success
!
use menu_f90
use gino_f90
      type (GPIXFRAME) frame
      type (GDIM) paper,newpaper
      character*(*) file
      integer,dimension(:),allocatable :: pix
      integer widget,istat,category,ndt,ndc,ip,metatype,xmax,ymax
      integer npix,l,iprop,nxpix,nypix
!     eps properties
      integer           :: properties(4) = (/1,4,3,0/)
!
!     first check widget is a graphics frame
!
    call gmActivateGraphicsFrame(widget )

      istat=0
      call gmEnqWidgetProp(widget,gmWidgetCategory=category)
      if(category.ne.GGRAPHICS .and. category.ne.GOPENGL) then
         call gmDisplayErrorBox('Save Graphics Frame', &
                          'Widget not a graphics frame')
         frasav=istat
         return
      end if
!
!     device must be in true colour mode for gGetPixelArea to work
!
      call gEnqColourInfo(ndc,ndt)
      if(ndt.lt.4) then
!         call gmDisplayErrorBox('Save Graphics Frame','Device not in true colour mode')
         write(6,*) 'Save Graphics Frame Device not in true colour mode'
         frasav=istat
         return
      end if
!
!     get frame dimensions
!
      call gmEnqGraphicsPixelBounds(widget,frame)
      xmax=frame%xmax
      ymax=frame%ymax
!
!     Allocate area for pixels
      npix=xmax*ymax
      allocate(pix(1:npix),stat=istat)
      if(istat.ne.0) then
         call gmDisplayErrorBox('Save Graphics Frame','Failed to allocate memory for pixel area')
         frasav=istat
         return
      end if
!
!     capture graphics frame in image
!
      call gGetPixelArea(0,0,xmax,ymax,1,1,xmax,ymax,pix)
!
!     parse file name to get extension (default = bmp)
!
      metatype=0
      l=gTrueLen(file)
     ! print *, 'Length=',l
      if(index(file,'.BMP').gt.0 .or. index(file,'.bmp').gt.0)&
                        metatype=1
      if(index(file,'.JPG').gt.0 .or. index(file,'.jpg').gt.0)&
                        metatype=2
      if(index(file,'.PNG').gt.0 .or. index(file,'.png').gt.0)&
                        metatype=3
      if(index(file,'.EPS').gt.0 .or. index(file,'.eps').gt.0)&
                        metatype=4
      if(index(file,'.WMF').gt.0 .or. index(file,'.wmf').gt.0)&
                        metatype=5
!
      if(metatype.eq.0) then
         if(index(file,'.').gt.0) then
            call gmDisplayErrorBox('Save Graphics Frame', &
                            'Unknown image file type')
            frasav=istat
            return
         end if
         file(l+1:l+4)='.bmp'
         l=l+4
         metatype=1
      end if
!
!     suspend guiwn and nominate image file device
!
      call gSaveGinoState() 
      call gSuspendDevice

      if(metatype.eq.1) call gBmp
      if(metatype.eq.2) call gJpeg
      if(metatype.eq.3) call gPng
      if(metatype.eq.4) call gEpsexp(1500.0,1500.0,0.0,0.0,4,properties)
      if(metatype.eq.5) call gWmfp(72,0,0,xmax,ymax)
      call gSetDeviceFilename(file(1:l),0)
      call gSetColourInfo(ndc,ndt)
!
!     set output image size to match graphics frame
!
      call gEnqDrawingLimits(paper,ip)
      call gEnqPixelResolution(nxpix,nypix)
      newpaper%xpap=float(xmax)/float(nxpix)*paper%xpap
      newpaper%ypap=float(ymax)/float(nypix)*paper%ypap
      call gSetDrawingLimits(newpaper,ip)
!
!     draw image to image file
!
      call gDrawPixelArea(0,0,xmax,ymax,1,1,xmax,ymax,pix)
!
!     deallocate pixel array
!
      deallocate(pix)
!
!     close image file and go back to gGuiwin/oglwin
!
      call gCloseDevice
      if(category.eq.GGRAPHICS) call gGuiwin
      if(category.eq.GOPENGL) call gOglwin
      call gRestoreGinoState(GMAPPED) 
      istat=1
      call gSetColourInfo(ndc,ndt)
!
      frasav=istat
      return
      end function frasav


end module etienne_gino
